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Abstract 

A new  concept  for  a retrodirective  target- 

✓ 

seeking  sonar  system  is  described.  The  system  is 
based  on  time  reversal.  This  is  achieved  by 
storing  the  target  return  in  a serial  access 
memory  and  reading  out  the  stored  information 
in  reverse  order  for  retransmission. 

The  theory  of  the  device  is  given  for  noise- 
free  operations  and  stationary  array  and  targets. 
The  motion  sensitivity  of  the  system  is  discussed 
and  the  effect  of  white  noise  and  reverberation 
noise  on  the  system  considered.  A computer 
simulation  of  the  system  in  a reverberation 
environment  is  described  and  some  results  of 
the  simulation  given. 
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I.  Introduction 

The  concept  of  retrodirective  and  target  seeking 
radar  or  sonar  systems  is  not  new.  In  1940  L.  C. 

Van  Atta  received  a patent‘d  for  a retrodirective 
phased  array.  It  functions  like  a corner  reflector 
sending  impinging  plane  waves  back  into  themselves. 
This  retrodirectivity  is  achieved  in  both  the  corner 
reflector  and  the  Van  Atta  array  by  inverting  the 
space  coordinates  of  the  wave.  In  the  Van  Atta  array 
the  space  inversion  is  done  by  connecting  the  receive 
array  elements  located  at  (x  , y ) to  the  corre- 
sponding  space  inverted  elements  (-xn,  -yn)  of  a 
transmit  array. 

If  one  puts  amplifiers  between  the  receive  and 
transmit  arrays  the  Van  Atta  array  can  be  made  target 
seeking.  A self-seeking  communication  system  based 
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on  this  concept  has  been  described  in  the  literature. 

We  shall  in  the  following  describe  a system  which 
is  based  on  time  inversion  instead  of  space  inversion. 
The  time  inversion  system  has  a number  of  desirable 
characteristics  which  cannot  be  achieved  with  a Van 
Atta  array.  These  characteristics  include:  self- 
focusing,  partial  motion  compensation,  target 
position  independent  echo  return  time  and  a tendency 
to  prevent  echoes  from  spreading  out  in  time. 

We  shall  discuss  in  this  report  possible 
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implementations  of  such  a system,  its  beamforming 
characteristics  in  a noise-free  environment  including 


Channel 

Processor 
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motion  sensitivity  and  configuration  stability,  the 
effect  of  signal  independent  white  noise  and  the 
system  performance  in  a reverberation  environment. 

II.  The  System 

The  concept  of  a time  inversion  retrodirective 
sonar  (or  radar)  is  straightforward;  it  is  based  on 
the  fact  that  an  outgoing  spherical  wave  changes 
under  time  inversion  into  an  incoming  spherical  wave 
converging  toward  its  center.  Therefore,  if  an 
array  of  transducers  receives  and  stores  over  a 
period  of  time  a signal  or  echo  from  a point  source 
and  retransmits  the  stored  signal  in  reverse  order, 
that  is,  inverted  in  time,  a beam  results  which  is 
focussed  on  the  original  source  position.  If  more 
than  one  source  is  present  in  the  field,  retro- 
directive  beams  focussed  on  each  one  of  the  sources 
are  generated.  This  characteristic  makes  the  time 
inversion  system  potentially  useful  as  a target 
seeking  and  tracking  device.  We  describe  in  this 
section  the  basic  components  of  such  a system  and 
give  a brief  summary  of  its  operational  character- 
istics. 

Figure  1 shows  a schematic  of  the  system.  It 
consists  of  an  array  of  N transducers,  each  connected 
to  its  own  channel  processor.  The  processor-trans- 
ducer combination  is  self-contained  and  the  N units 
can  be  distributed  arbitrarily  in  space,  connected 
only  by  wire  or  radio  link  to  a common  clock  and 


5 


display.  In  one  possible  implementation  the  central 
element  of  the  array  is  connected  to  a search  pulse 
generator.  The  target  seeking  process  in  this  case 
is  initiated  by  the  emission  of  a search  pulse  by 
this  central  element.  We  shall  use  this  model  for 
our  discussion  of  the  system  characteristics. 

Figure  2 shows  one  transducer  with  its  channel 
processor.  A clock  operates  the  transmit-receive 
gate,  samples  the  data  and  controls  their  transport 
in  and  out  of  a serial  access  memory.  In  the  receive 
mode  the  signal  is  amplified,  sampled  and  clocked 
into  a serial  access  memory.  The  "store-read  gate" 
routes  the  data  samples  in  and  out  of  the  memory  and 
supplies  appropriately  coded  information  to  a display 
unit  which  displays  the  content  of  the  entire  memory 
bank. 

The  system  performance  depends  critically  on  the 
characteristics  of  the  memory  device;?.  The  number  N 
of  stored  samples  which  are  necessary  to  adequately 
describe  a signal  of  duration  T and  bandwidth  B is 
according  to  the  Nyquist  sampling  theorem: 

N = 2BT  (1) 

The  signal  duration  T is  directly  related  to  the 
range  R over  which  targets  are  being  acquired  by  the 
system: 

T = 2R/c  (2) 

where  c is  the  signal  propagation  velocity.  The 
bandwidth  B is  similarly  related  to  the  range 
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resolution  AR: 
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B = c/2AR  (3) 

The  maximum  range  R and  the  resolution  AR  of  a time 
inversion  system  are  limited  by  the  retention  time 
AT  of  the  memory  and  the  maximum  clock  rate  p at 
which  data  samples  can  be  moved. 

Because  of  the  time  inversion  the  first  signal 
sample  clocked  into  the  memory  is  the  last  one  to  be 
clocked  out  and,  therefore,  stays  in  the  memory  for  a 
period  of  time  equal  to  twice  the  signal  return  time 
T.  We  have,  therefore,  a limiting  condition  for  the 
retention  time  AT: 

AT  > 2T  (4) 

and  because  of  Equation  (2)  a limiting  condition  for 
the  range  R: 

R < cAT/4  (5) 

The  minimum  clock  rate  necessary  to  store  an  incoming 
signal  of  bandwidth  B is: 

p . = 2B  (6) 

min 

which  gives  together  with  equation  (3)  the  limiting 
condition  for  the  range  resolution  AR: 

AR  = c/p  (7) 

If  we  consider  a typical  high  resolution  sonar 
with  a search  range  R of  500  meters  and  a range  reso- 
lution AR  of  .1  meter,  we  find  from  the  foregoing 
considerations  the  following  requirements:  A storage 
capacity  of  N = 10  bins,  a minimum  retention  time  of 
T = 800  msec  and  a clockrate  of  p = 25  kHz.  These 
requirements  have  to  be  met  together  with  a high  dynamic 
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range  and  reversibility  of  the  dataflow  in  a serial 
access  memory.  Presently  available  commercial  charge 
coupled  devices  (CCD)  or  similar  devices  which  satisfy 
the  reversibility  requirement  almost  meet  these 
criteria.  The  most  critical  parameter  is  the  reten- 
tion time  AT.  Since  the  dynamic  range  of  CCD's 
deteriorates  with  retention  time  and  ambient  temper- 
ature, the  above  time  requirement  can  be  met  with 
presently  available  CCD's  only  if  they  are  cooled. 

The  required  number  N of  storage  bins  is  high  but 
feasible,  and  the  required  clockrate  orders  of 
magnitude  below  device  capability. 

If  we  consider  a search  radar  as  a potential 
application  and  assume  a range  of  50  km  and  a range 
resolution  of  10  meters  we  find  that  the  same  memory 
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capacity  of  N = 10  storage  bins  is  required,  but  a 
much  smaller  retention  time  of  T = .8  msec  and  a 
minimum  clockrate  of  p > 30  MHz.  In  this  case  the 
clockrate  is  the  critical  parameter  but  still  within 
the  capability  of  CCD's.  The  expected  technical  and 
economic  feasibility  of  an  adequate  memory  bank  for 
a time  inversion  system  in  the  near  future  makes  it  a 
timely  task  to  analyze  its  characteristics. 

III.  The  Array  Target  Interaction 

In  this  section  we  derive  the  algebra  which 


describes  the  interaction  of  the  time  inversion 
system  with  a field  of  targets. 
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We  assume  an  array  of  K transducers  and  a field 

of  J taryets.  The  positions  of  the  transducers  and 

targets  are  described  by  position  vectors  in  a 

cartesian  coordinate  system  (x,y,z)  which  has  its 

origin  in  one  of  the  array  elements,  called  the 

reference  element.  The  reference  element  is  the 

element  closest  to  the  center  of  gravity  of  the  array. 

It  is  assumed  that  the  radius  r . of  the  smallest 

min 

sphere  which  contains  all  array  elements  is  small 

compared  to  the  minimum  distance  R . from  which 

min 

targets  are  acquired  by  the  system.  The  maximum 

distance  from  which  targets  are  acquired  is  R 

^ max 

The  range  R of  operation  of  the  system  is  R-R  -R  . . 

c max  min 

We  limit  the  discussion  in  this  section  to  a 
linear  system,  that  is,  we  assume  that  the  n+lst 
retransmitted  signal  is  a linear  function  of  the  nth 
target  response  which  is  in  turn  a linear  function 
of  the  nth  retransmitted  signal.  It  is  the  purpose 
of  this  section  to  derive  the  linear  operators  which 
relate  array  output  and  target  response  to  each  other. 

In  order  to  avoid  unnecessary  complication  we 
assume  that  the  linear  dimension  'a'  of  the  individual 
transducers  are  of  the  order  of  the  mean  wavelength 
of  the  (bandlimited)  search  function.  The  farfield 
pressure  wave  generated  by  the  kth  element  during  the 


nth  retransmission  is,  therefore,  within  a cone  angle 

nrk . 

, a snh.->rioal  wave* 

a 

n V f 

characteristic  function  f (t)  of  this  spherical  wave 


of  0 = sin_1(^-)  a spherical  wave  — - — ..r/-c)  . The 
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has  the  dimension  [Newton  per  meter]  and  is  related 
to  the  voltage  nV^  across  the  kth  transmitting  trans- 
ducer by 

nfk  - »V  (8) 

with  the  coefficient  a given  by 

a = { (4a2npc)/(TrXI)  }1/2  (9) 

where  n is  the  efficiency  of  the  electrical  to 
mechanical  energy  conversion,  p the  density  of  the 
transmitting  medium,  and  I the  electrical  impedance 
of  the  transducer. 

We  describe  the  total  array  output  for  the  nth 
retransmission  as  the  ordered  set  nF(t)  of  character- 
istic functions  nfk(t): 

nF(t)  = (nfk(t)}  k = 0, K— 1 (10) 

We  assume  further  that  the  targets  are  "point 
scatterers"  in  the  limited  sense  that  they  are  small 
compared  to  the  lateral  resolution  of  the  array.  This 
implies  that  the  backscattered  wave  is  again  adequately 
described  by  a spherical  wave  nhj (t-r/c)/r  within  the 
small  cone  which  the  array  subtends  at  the  target.  The 
target  response  function  nH(t)  is  then  again  the 
ordered  set  of  the  generating  functions  nhj (t)  of  these 
spherical  waves: 

nH(t)  = {"hjtt)}  j = 1 J . (11) 

The  limitation  of  our  consideration  to  point  targets 
does  not  limit  the  generality  of  our  derivation  since 
one  can  build  up  arbitrarily  complex  targets  by  an 
ensemble  of  suitably  chosen  point  targets  (Huyghens 
Principle) . 
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It  follows  from  the  foregoing  definitions  that 


the  pressure  npj (t)  at  the  jth  target  location  due 

to  the  nth  retransmission  nF(t)  is 
n 


" k^nfk<t"V 


(12) 


where  r..  is  the  distance  between  the  jth  target  and 

the  kth  transducer  and  tjk  = r^/c  the  corresponding 

propagation  time.  It  is  sufficient  to  approximate 

in  the  slow  varying  spreading  term  1/r^  by 

r.  = r.  , the  distance  of  the  jth  target  from  the 
J D o' 

array  center. 

We  introduce  the  time  delay  operator  D(t) 
defined  by 

D(x)f(t)  = f(t-x)  (13) 

and  write  equation  (12)  as: 


np.(t)  = |-£  D(tjk)nfk(t) 
J K 


(14) 


The  jth  target  response  is  a linear  function  of 


(15) 


the  pressure  field  at  all  target  locations 

nh.  (t)  = I S.  V,  (t) 

J j I J J J 

The  ordered  set  {Sjj,}  is  the  scatter  matrix  of  the 
target  distribution.  If  we  neglect  multiple  scattering 
the  target  response  depends  only  on  the  primary  pres- 
sure field  nPj (t)  at  its  own  location.  The  sum  in 
equation  (15) , therefore,  reduces  to: 

nhj  (t)  = Sjtt)  Vj(t)  = Ojijtt)  * nPj(t)'  (IS) 

where  0j  is  the  backscatter  cross  section  and  ij (t) 
the  normalized  impulse  response  of  the  jth  target. 


MM 
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The  star  in  eq.  (16)  indicates  convolution.  Intro- 
ducing equation  (14)  into  equation  (16)  we  obtain: 

"V«  = Vrj  * l D(tjk>n£k(t> " 


Vr3  l D(Vi:i(t)  * £k<£> 


To  go  to  the  right  side  of  eq.  (17)  we  made  use  of 
the  commutability  of  the  convolution  operation  with 
the  time  delay  operator.  This  follows  directly  from 
the  definition  of  the  convolution  operation: 


g (t)  * f(t) 


g(t) f (t-x)dr 


We  have  shown  above  that  in  the  neighborhood  of 
the  array  the  waves  scattered  back  from  the  targets 
can  be  approximated  by  zero  order  spherical  waves. 

We  obtain,  therefore,  for  the  pressure  npk(t)  at  the 
location  of  the  kth  transducer: 

"pk(t)  = ? FT  "hj (19) 

The  voltage  nv^(t)  generated  by  the  kth  receiving 
transducer  due  to  the  pressure  npk(t)  is: 

nvk(t)  = enpk(t)  (20) 

where  B is  the  transducer  constant. 

In  order  to  determine  the  next  retransmitted  signal 
n+1fk(t)  we  have  to  keep  track  of  the  appropriate 
receive  and  transmission  times,  which  are  shown 
schematically  in  Fig.  3.  We  label  the  time  at  which 
the  nth  receive  period  is  activated  as  nta  and  the  time 
at  which  it  ends  as  nt^.  During  the  receive  time 
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nta  < t < nt^  the  transducers  generate  the  voltages 
nvk(t)  which  are  amplified,  filtered,  sampled  and 
stored.  At  t = n+^ta  = nt^  the  n+lst  transmit  mode 
is  activated.  The  stored  signals  are  read  out  in 
reverse  order,  filtered,  amplified  by  a large  factor  y 
(power  amplification)  and  fed  to  their  respective 
transducers.  The  transmit  voltage  xVk  across  the 
kth  transducer  becomes : 

n+1Vk(t)  = Y nvk{nte-(t-ntB) } (21) 

and  the  retransmitted  signal  n+'*'fk(t)  according  to 
eq.  (8) 

n+lfk(t)  = all+lvk(t)  (22) 

The  retransmission  ends  when  all  of  the  stored  signal 


is  retransmitted  that  is  at  time: 

t = n+1tQ  = n+1t  + (nt^  - nta)  = 2nt^  - nta  (23) 

p ct 

Combining  equations  (19)  through  (22)  we  obtain: 

n+1fk(t)  = D(2nte)S  E>(tjk)nhj  (-t)  (24) 

where 

K = agy  (24a) 

Equations  (17)  and  (24)  are  sufficient  to  relate  the 
n+lst  retransmitted  signal  back  to  the  original  search 
pulse. 

We  assume  for  the  following  analysis  that  the 
initial  pulse  is  transmitted  by  a single  projector  at 


the  location  of  the  array-element  k = 0.  The  transmit 
gate  opens  at  t = 0tQ  and  closes  at  °tg.  In  order  to 
accommodate  all  later  even  retransmissions  we  let: 
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where  r . is  the  radius  of  the  smallest  sphere 
min 

capable  of  enclosing  the  array  and  t the  pulse 

duration.  The  transmit  time  given  by  eq.  (25)  is 

the  minimum  time  required  to  allow  beamforming  in 

any  desired  direction  with  all  array  elements 

participating.  The  projector  which  is  in  the 

center  of  the  array  in  our  model  transmits  during 

the  time  °t  + r . /c  - t/2  < t < °t  + r_.  /c  + x/2. 
a min  a min 

The  target  response  to  the  initial  search  pulse 
°f(t)  is  according  to  equation  (16) 

ci 

“hjtt)  = D(tjo)ij  * f(t)  ' (26) 

In  order  to  receive  at  all  transducer  locations  the 

backscattered  signal  from  all  potential  targets  in 

the  acquisition  range  R,  it  is  necessary  for  the 

receive  mode  to  begin  at: 

t = °ta  = 2(R  . - r . )/c  + °t  (27a) 

' min  min"  a 


and  to  run  until: 

t = °t“  = 2 (R  . - r • )/c  + °t  = 

min  min"  a 

-(R  + 2r  . ) + 2t  + °t 

c'  max  min  a 


(27b) 


All  potential  target  returns  from  within  the  accep- 
tance range  R of  the  system  are  now  received;  the 
array  can,  therefore,  begin  retransmission  imme- 
diately. Labeling  the  beginning  of  the  first 
1 18 

retransmission  ta  we  have:  tft  = #t  . The  trans- 

mission cycle  runs  until  stored  information  is 
retransmitted,  that  is,  until 

t - % - 20te  - % - 

f<2!W  - ‘Win  + 5rmi„>  + + 


(28) 
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The  first  retransmitted  signal  f^(t)  is 
obtained  by  introducing  equation  (26)  into 
equation  (24) : 

1 6 ^ <o-D(t--t-.  ) 

•Lfv(t)  = D(2  °tP)  E —3 L-JJL-  i (_t)*f(-t)  (29) 

k 1 j-1  (r.)*  3 

The  target  response  1hj (t)  to  that  excitation 
follows  if  we  introduce  equation  (29)  into  equation 
(17) . After  some  rearrangement  of  the  resulting 


double  sum  we  obtain: 

6, 


hj(t)  = — r 


<a-D(2  °t'J) 


E E 
k j'  (r j , ) 


2 D<tjk”tj'k"tj') 


(30) 


ij  (t)  *i  . , (-t)  x f (t) 


We  note  that  the  term  with  j ' = j in  the  sum  over  j ' 
in  eq.  (30)  corresponds  to  the  excitation  of  the  jth 
target  by  its  own  retrodirected  beam.  The  contribu- 
tions from  all  k transducers  become  equal  and  add 
coherently.  The  terms  with  j ' ^ j are  the  side  lobe 
contributions  at  the  jth  target  from  beams  focussed 
at  the  other  targets.  To  set  this  in  evidence  we 
rewrite  equation  (30)  to  read: 

1 ft  KKai2 

Ah.(t)=D(2  °tp)«{ D(-t.)a. (t)*f (-t)  + 

3 (r  j ) 3 3 3 


D(t..-t.  ,.-t.  ,)K-0., 

E E JK  V ' J 1-3-  • c.  .,(t)*f(-t)} 


(31) 


k j Vj 


r3  rj' 


33 


where  a^ (t)  is  the  autocorrelation  function  of  the 
target  impulse  response: 


aj(t)  = ij  (t)*ij(-t) 


(32) 
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and  Cjj,(t)  the  cross  correlation  between  the  impulse 
responses  of  the  jth  and  j'th  target: 


cjjt  (t)  = ij(t)*ij, (-t) 


(33) 


The  side  lobe  contributions  to  xh(t)  given  by 
the  double  sum  in  equation  (31)  are  in  general  very 
small  compared  to  the  focussed  beam  contributions 
given  by  the  first  term  in  eq.  (31) . This  is  true 
firstly  because  the  side  lobe  level  of  pulsed  systems 
is  already  low  and  secondly  because  the  cross 
correlation  functions  Cjj,  of  the  impulse  responses 
of  different  targets  can  be  expected  to  be  small 
compared  to  their  autocorrelation  functions.  However, 
another  very  effective  side  lobe  reduction  specific 
to  the  present  system  is  possible  through  appropriate 
timing  of  the  receive  period  during  which  the  back- 
scattered  signal  ^hj (t)  is  accepted. 

The  pressure  field  ^p^ft)  at  the  array  due  to 

1hj (t)  follows  from  eqs.  (19)  and  (31): 

2 

1 - ft  KO  . K 

pk(t)  = D (2  °ty){£  — ^-4  DU^-tjJa.  (t)*f(-t)  + 


j <rj) 


(34) 


+ Z £ 

k jVj  ( 


^7-72  D(tjk,+tjk"tj,k,-tj,)Cjj'(t^£(“t)} 


y ,-3 

Equation  (31)  shows  that  the  main  lobe  contribu- 
tions from  all  targets  arrive  at  the  reference  element 
(k=0)  at  the  same  time  t = 2 °t®.  If  the  targets  are 
point  targets  the  autocorrelation  functions  a^ (t)  are 
delta  functions  and  the  backscattered  field  seen  by 
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i 


the  array  elements  is  proportional  to  the  search 
pulse  (inverted  in  time).  Due  to  the  finite  ex- 
tension of  the  array,  the  first  signal  at  any  array 

element  arrives  at 

.3 


t > 2 °t  - r . /c 
- u min/c 


6 


and  has  passed  the  array  at  t < 2 °t  + r . /c  + t. 

- min 


If  one,  therefore,  sets  the  receive  gate  to  open  at: 


1ta  = 2 °tB  - r . /c  = - R 
min'  c m, 


+ — r + 4t  + °”t. 
ax  c rain  a 


(35) 


and  close  at: 


t = 1t^  = 2 °tP  = 2 °te  + + 

c 


t = — R + g/c 
c max 


r . + 5t  + °t 

min  a 


all  main  lobe  signal  returns  will  be  received.  For 
extended  targets  for  which  the  autocorrelation 


functions  a^ (t)  are  not  delta  functions,  but  by 


definition  symmetrical  functions  with  a maximum  at 
t = 0,  the  same  receive  gate  setting  still  assures 
the  reception  of  a major  portion  of  the  main  lobe 
return  signals  from  all  targets. 


1.3 


The  short  duration  of  the  receive  mode  ( t - 


1t“  = T + 


2r 


min 


) excludes  most  side  lobe  returns. 


They  are  accepted  only  from  target  pairs  with  a range 


difference  rjji  of  the  order  of  the  linear  array 


dimension: 


r . . , < 4r  . + tc 

33 1 - min 


(36) 


Equation  (27a, b)  and  (28)  assure  that  the 


beginning  of  the  receive  mode  ^ta  is  compatible  with 


the  end  of  the  previous  transmit  node. 
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After  all  main  lobe  returns  are  received,  the 

system  is  ready  for  the  next  transmission  period  which 
2 18 

starts  at  t = t = t.  The  retransmission  lasts  till 
all  stored  signals  are  retransmitted,  that  is,  until: 
t = 2tB  = 2ta  + 1te  - 1ta  = 21t&  - ^-t01. 


The  retransmitted  signal  ^f^(t)  now  follows  from 
eqs.  (31)  and  (24): 

2 2 

2 1 * K 

fk(t)  = D(2±tft)  Z V D(t.-t,.  )a.(t)*f(t)  + 

k 3 j (r,)4  3 3 (37] 

J side  lobe  terms. 

The  side  lobe  terms  are  small  and  almost 

always  gated  out.  We  shall,  therefore,  neglect 

them  from  here  on.  Note  that  the  autocorrelation 

function  a(t)  is  not  affected  by  time  reversal 

since  it  is  by  definition  a symmetric  function. 

The  time  delay  factor  D(21t^)  in  front  of  the  sum 

can  be  deleted  if  we  introduce  a new  time  reference: 

t'  = t - 21tg  (38) 

and  write  the  running  time  again  as  t instead  of  t'. 

Equation  (37)  implies  that  the  total  array  output 

for  the  second  retransmission  represents  J beams 

focused  back  onto  the  J targets.  This  insonification 

2 

generates  the  target  response  hj (t)  which  we  obtain 
by  substituting  eq.  (37)  into  eq.  (17)  again 
neglecting  all  side  lobe  terms: 


h.(t)  = 


9,a.D(t.) 


aj (t) *ij (t)*f (t) 


with 


kKct  . 

gj  = — i 


(39a) 
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If  one  compares  (t)  given  by  eq.  (39)  with  the 
target  response  °hj (t)  due  to  the  initial  search 
pulse  which  is  given  in  eq.  (26)  one  obtains: 

2hj(t)  = gj2aj(t)*°hj(t)  (40) 

The  autocorrelation  function  a^ (t)  of  the  target 
impulse  response  is,  up  to  a scale  factor,  a 6 -function 
for  point  targets  or  a narrow  pulse  like  function  for 
complex  extended  targets.  Therefore,  eq.  (40)  implies 

that  the  second  hand,  therefore,  all  even)  iterations 

> 

of  the  target  response  are  more  or  less  perfect 
replicas  of  the  target  response  to  the  initial 
search  pulse: 

2nhj (t)  = gj2  • aj(t)*2(n_1)h.t  = g2najn(t)*°hj (t)  (41) 

where  a^n  is  the  n-th  iteration  of  the  a^ (t) : 

ajn (t)  = aj(t)*aj(t)  *a j (t)  (42) 

Since  the  n-th  target  output  nfk(t)  is  a linear 

St 

function  of  the  (n-1)  target  response,  the  quasi- 
periodicity of  nhj  carries  over  to  nf (t)  : 

2n  gi2nD(_tki) 

fk(t)  = ? ajn(t)*f(t)  (43) 


and 


2n+l 


2n+l£k“>  - £ 


D(tkj-tj)oj(-t)*ajn(t)*f  (-t)  (44) 


If  the  round-trip  gain  factor  g^  equation  (43)  is 
greater  than  unity  the  array  output  increases  with 
every  iteration  (equation  43,44).  If  the  round-trip 
gain  factor  for  a given  target  is  smaller  than  unity 
then  the  target  insonification  dies  out.  Since  g is 
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proportional  to  the  scatter  cross  section  and 
inversely  proportional  to  the  target  distance,  a 
free-running  system  singles  out  the  strongest  and 
closest  targets.  If  the  output  signal  amplitude  is 
limited,  the  system  tends  to  automatically  equalize 
the  beam  strength  for  each  target  beam  provided  g^ 
is  greater  than  unity. 

Equations  (43)  and  (44)  show  two  significant 
differences  between  the  even  and  odd  retransmission 
functions.  The  even  retransmission  represents  j 
simultaneously  transmitted  beams  focused  on  the  j 
targets.  The  waveform  on  each  beam  is  practically 
identical,  consisting  of  a more  or  less  true  replica 
of  the  waveform  of  the  search  pulse.  In  contrast  the 
odd  retransmissions  represent  j beams  which  are 
spaced  out  in  time  and  which  have  target  specific 
waveforms.  The  waveforms  are  essentially  the  initial 
pulse  convolved  with  the  respective  target  impulse 
response.  They  preserve,  therefore,  the  target 
signature.  The  emission  time  of  each  beam  is 
proportional  to  the  range  of  its  target.  This  is  a 
useful  feature  since  it  permits  one  to  select  one 
or  more  particular  targets  by  simply  gating  out  the 
corresponding  beams.  This  gating  has  to  be  done  only 
once;  from  there  on  the  system  continues  to  track  only 
the  selected  targets. 

We  arrived  at  the  relatively  simple  results  of 
eqs.  (43)  and  (44)  by  neglecting  all  side  lobe 

| 

A 
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contributions.  To  verify  the  validity  of  this 
approximation,  a computer  model  of  the  system 
including  all  side  lobe  contributions  was  developed. 
The  program  is  given  in  Appendix  I.  Some  results 
of  the  model  study  are  shown  in  Figs.  4 and  5.  The 
following  model  was  used  to  obtain  this  result: 

The  transmit  receive  unit  of  the  system  was  a linear 
array  of  513  equally  spaced  elements  with  d = X/2. 
The  center  element  was  at  k = 0.  It  emitted  the 
initial  search  pulse  °f(t).  The  search  pulse  was 
chosen  to  be  a single-frequency  pulse  with  a 
gaussian  envelope: 


_/Vt) 

°f(t)  = e“l2lTVt  e 5 

The  system  operated  against  a field  of  5 point 
targets  of  equal  cross  section.  The  location  of  the 
targets  is  shown  in  Table  1. 


Table  1 


Target 

z/103X 

1 

-6 

20 

2 

-2 

20 

3 

0 

20 

4 

+2 

20 

5 

+6 

20 

Target  Locations 


J 
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The  target  positions  were  chosen  symmetrically  around 

the  z axis  to  assure  that  side  lobe  contributions  were 

. 

received  during  the  short  receive  period  At  = 

(^"tg  - ^ta)  . Figure  6 shows  the  envelope  of  the 
output  voltage  ^v^.(t)  eq.  (24)  of  the  receiving 
array  elements  during  At  for  the  element  k = 2.  •'»  <n 
a logarithmic  scale  which  puts  the  side-lobe 
contributions  in  evidence.  Note  that  the  side-lobe 
contributions  are  two  orders  of  magnitude  smaller 
than  the  main  returns.  Figure  7 shows  the  returns 
of  elements  256,  192,  128,  69,  and  0 on  a linear 
scale.  The  side-lobe  effects  are  no  longer  visible. 

The  beamforming  delays  between  channels  can  be 
. clearly  seen. 

IV.  Motion  Sensitivity 

The  time  inversion  systems  described  above 
focus  the  energy  of  each  successive  ping  back  to  the 
area  where  the  previous  ping  encountered  a target. 

The  effect  of  this  self-focusing  is  to  enhance  the 
target  insonif ication  in  successive  pings  over  that 
obtainable  with  a non-directional  system.  This 
objective,  however,  can  be  reached  omly  if  the  target 
has  not  moved  out  of  the  focal  area  diuring  the  time 
between  pings.  Since  the  linear  dimension  of  the 
focal  area  is  inversely  proportional  to  the  aperture 
of  the  array,  this  motion  sensitivity  imposes  a limit 
to  the  useful  aperture  of  a system  operating  against 
moving  targets.  Since  the  lateral  estent  of  the  focal 

^ 1 J 
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area  is  smaller  than  the  depth  of  focus,  the  critical 
target  motion  is  the  motion  perpendicular  to  the 
line  connecting  the  target  with  the  array.  In  order 
to  find  the  maximum  aperture  size  which  can  be  used 
effectively  against  moving  targets,  we  need  to 
consider  only  this  worst  case.  If  the  time  between 
pings  is  At  and  the  target  velocity  is  v,  the  distance 
travelled  between  pings  by  the  target  is 

d = At*v  (45) 

This  distance  has  to  be  compared  with  the  halfwidth 
of  the  focal  area  F: 


where  X is  the  wavelength  of  the  center  frequency  of 
the  search  pulse  and  D the  effective  diameter  of  the 
array  aperture.  Since  by  design  the  target  starts  out 
in  the  center  of  the  focal  area,  the  travelled  distance 
d between  pings  must  be  not  greater  than  1/2  F.  The 
limiting  condition  for  the  array  aperture,  therefore, 
becomes : 


XR. 

d i 2d1 


(47) 


or 


XR. 


D i Dc  - 231  - 


XRi 


2At*  v 


(48) 


where  Dc  is  the  critical  aperture.  He  have  seen  above 
that  the  time  between  pings  in  our  system  depends  on 
whether  the  last  ping  was  due  to  an  odd  or  even  re- 
transmission. Referring  to  Fig.  (3)  and  eq.  (27b)  we  find: 


..  2n.f5  2n.a  , 

c’Atodd  = t - t + 


2Rm 


ax 


- 2R 


J = 4Rn. 


ax  - 2Rj 


(49) 


I 
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Equation  (49)  shows  that  for  odd  retransmission 
the  elapsed  time  between  the  present  and  the  previous 
ping  is  largest  for  targets  close  to  the  array. 

This  is  due  to  the  fact  that  the  returns  from  nearby 
targets  spend  extra  time  in  the  system  memory.  The 
elapsed  time  for  even  retransmissions  is  simply: 

°'Ateven  * 2"+lt6  ' 2"+lt“  ' 2Rmax  + 2Rj  = 2Rj  <50> 
which  is  for  all  targets  smaller  than  AtQdd  of  eq.  (49) . 

Since  the  system  has  to  accommodate  the  largest  elapsed 

time  (At,,)  : 

odd  max 

C*  (At  ..)  = 4R  -2R  . = 2 (2R+R  . ) 

odd  max  max  mm  mm 

we  obtain  for  the  critical  aperture  D of  eq.  (48) : 


D ■ l..C..Vn 

C 4v ( 2R+R 


mm 


(51) 


The  severity  of  the  aperture  limitation  given  by 
eq.  (51)  depends  on  the  factor  c/v.  In  the  radar 
case  where  the  factor  c/v  for  target  speeds  of  interest 
is  very  large  the  critical  aperture  D is  much  larger 
than  the  aperture  size  required  for  a practical  system. 

For  target  velocities  of  Mach  10  the  factor  c/v  is  still 

5 4 

about  10  and  D , therefore,  of  the  order  of  10  X. 
c 

The  situation  is  quite  different  in  the  sonar  case. 
Here  c/v  for  targets  with  a velocity  of  50  km/h  is 
only  about  10  and  the  critical  aperture,  therefore, 
of  the  order  of  a few  wavelengths,  which  is  too  small 
for  a practical  search  and  tracking  sonar.  We  shall 
discuss  in  the  next  section  a simple  system  modification 
which  partially  overcomes  this  difficulty. 
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System  Modification  for  Tracking  Moving  Targets 
The  critical  aperture  Dc  equation  (51)  is 
required  to  assure  proper  insonif ication  of  all 
targets  during  all  retransmissions.  The  limit  is 
imposed  by  the  odd  retransmissions  and  within  the 
odd  retransmissions  by  the  targets  closest  to  the 
array. 


For  even  retransmissions  the  permissible 
aperture  5 is  independent  of  the  target  range  and 
equal  to: 

5 = d/(c*Ateven)  = Xc/2v 


(52) 


The  aperture  required  to  properly  insonify  the  most 
distant  target  during  odd  transmissions  is  equal  to  D 


since  for  R ■ = 

3 max 


the  time  lapse  At0^  becomes  equal 


to  At 


even 


c*  At 


, , = 2 (2R  -R- ) = 2R ■ = c * At  for  R,  = R 

odd  max  j j even  j max 


(53) 


One  can,  therefore,  satisfy  the  insonif ication  require- 
ment for  distant  targets  by  starting  out  the  odd  re- 
transmit period  with  the  full  aperture  5 (which  is 

4R 

greater  than  D by  a factor  (2  + = ) and  reduce  the 

c min 

aperture  linearly  with  time  to  reach  Dc  at  the  end 
of  the  odd  retransmit  period.  The  full  aperture  5 is 
used  for  the  even  retransmits.  This  modification  can 
be  made  simply  by  adjusting  the  transmit  gates  for  the 
different  array  element  such  that  toward  the  end  of 
the  odd  transmit  period  only  a central  portion  of  the 
array  is  left  to  transmit. 
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The  net  effect  of  this  system  modification  is 
that  the  diameter  at  focus  for  all  beams  is  commen- 
surate with  the  maximum  distance  travelled  by  the 
generating  target. 

The  most  severe  loss  in  angular  resolution  and 
beam  intensity  occurs  for  targets  closest  to  the 
transmitting  array,  that  is,  for  the  case  where  it 
is  least  detrimental  since  for  the  closest  target 
the  aspect  angle  is  largest  and  the  spreading  loss 
of  the  beam  a minimum. 

4R 

In  spite  of  the  improvement  by  a factor  (2  + ^ ) 

min 

the  new  maximum  aperture  5 is  still  relatively  small  and 
represents  for  sonar  applications  a practical  system 
limitation.  For  target  velocities  up  to  50  km/h, 
for  example,  the  maximum  useful  aperture  is  still 
only  50  X. 

We  shall  show  in  the  following  section  that  this 
aperture  limitation  can  be  removed  completely  for  a 
system  which  operates  from  a moving  platform  against 
a field  of  stationary  targets. 

Moving  Array  Stationary  Targets 

We  consider  in  this  section  the  practically 
important  case  of  an  array  mounted  on  a moving  platform 
operating  against  a field  of  stationary  targets.  A 
mine  hunting  sonar  is  an  example  of  such  a system.  We 
shall  show  that  a minor  system  modification  removes  in 
this  case  all  aperture  limitation. 


The  aperture  limitation  discussed  above 


guaranteed  that  the  radius  at  focus  of  each  beam  was 
larger  than  the  maximum  distance  a target  could 
travel  between  pings.  This  is  necessary  because  in 
this  system  the  retrodirected  beam  is  focussed  on  the 
spot  where  the  last  ping  encountered  a target.  In 
the  moving  platform  situations,  one  can  use  the 
knowledge  of  the  platform  velocity  to  partially 
compensate  for  the  apparent  target  motion  by  advancing 
appropriately  the  focal  point  of  the  retrodirected  beam. 

We  consider  first  the  even  retransmissions. 

* 

Complete  time  reversal  would  require  that  the  platform 
moves  during  the  transmission  period  backwards  along 
. the  same  path  it  moved  during  the  receive  period. 

The  re-emitted  beam  would  then  focus  exactly  on  the 
target  location.  One  can  simulate  the  effect  of  the 
reverse  platform  motion  on  the  transmitted  wave  by 
slowing  down  by  a factor  (l-2v/c)  the  clock  rate 
which  controls  the  data  readout  during  retransmission. 
The  factor  (l-2v/c)  contains  twice  the  doppler  cor- 
rection v/c.  This  comes  about  since  we  want  to  simulate 
reverse  motion  on  a platform  which  continues  to  move 
forward.  The  retransmitted  beam  is  now  focussed  back 
onto  the  target  with  the  focal  point  displaced  only  by 
the  small  distance  A = (2nt^  - 2n-1ta)v,  which  the 
• platform  has  moved  between  the  beginning  of  the  recep- 

tion and  the  end  of  transmission  periods.  The  dis- 
placement A is  of  the  order  of  D»v/c,  that  is,  even  for 
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the  closest  targets  small  compared  to  the  focal 
XRmin  ' 


diameter 


In  other  words  we  can  for  the  even 


transmission  keep  the  target  in  focus  by  slowing  down 
the  transmission  clock  rate  by  the  factor  (l-2v/c) . 

The  same  argument  holds  for  odd  retransmissions  if 
the  target  is  at  the  end  of  the  range  Rj  = RJnax  . 

The  changing  aperture  approach  discussed  in  the 
previous  section,  therefore,  applies  also  to  the 
moving  platform.  No  restrictions  on  the  maximum 
aperture  are  now  necessary.  The  only  requirement 
is  to  make  the  clock  rate  for  all  retransmissions 
smaller  by  a factor  (l-2v/c)  than  that  for  receptions. 

Internal  Array  Motion  for  Non-rigid  Arrays 

If  the  array  consists  of  independent  units  as 
suggested  in  Section  2 (deployed,  for  example,  as 
sonobouys) , the  most  critical  motion  affecting  the 
operation  of  the  system  is  the  displacement  of  the 
array  elements  against  each  other  during  the  receive- 
transmit  times.  These  time  intervals  are,  as  we  have 
seen  earlier  in  Section  II,  orders  of  magnitude 
different  for  the  odd  and  even  receive-transmit  periods. 

At  = (2ntP  - 2n-1t  ) = 4r  . _/c  + x - 2D/c 
even  ' a min' 

4todd  ‘ <2n+lt6  - 2"V  ‘ 4R''C  <54) 

At-.  is  of  the  order  of  seconds;  At  is  of  the 
odd  even 

order  of  milliseconds.  The  permissible  relative 
displacement  of  array  elements  between  reception  and 
retransmission  is  a fraction  of  a wavelength  for 
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coherent  signal  addition  at  the  target.  This 

condition  is  not  generally  satisfied  for  the  long 

time  interval  AtQdd.  One  rather  has  to  expect  that 

for  higher  resolution  sonars  the  array  element 

position  is  a random  variable  with  a variance 

approaching  or  exceeding  a wavelength.  For  the 

short  time  interval  At  , however,  one  can  assume 

even 

that  coherence  between  received  and  transmitted  beam 
is  maintained. 

This  implies  that  for  odd  retransmissions  the 
intensities  of  the  various  array  element  contributions 
add  up  at  the  targets,  whereas  for  even  retransmission 
the  amplitudes  add  coherently  at  the  target  to  give 

v 

the  full  array  gain.  We  have  shown  earlier  that  the 
round-trip  gain  factor  g defined  in  equation  (40)  is 
a useful  parameter  to  describe  the  system  performance. 
The  round-trip  gain  factor  g&  for  the  coherent,  even 
retransmissions  is  equal  to  the  round-trip  gain  g 
given  in  eq.  (40) ; for  the  incoherent,  odd  retrans- 
mission the  round-trip  gain  factor  godd  becomes  due 
to  the  addition  of  intensities,  rather  than  ampli- 
tudes, approximately  equal  to: 

(55) 

The  average  gain  factor  g for  a full  cycle  of  even 
and  odd  retransmissions  for  the  non-rigid  array  is 
therefore: 

g2  - ge  • gQ  = g2/''*  (56) 


»odd  ■ 9A* 
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Since  the  gain  factor  for  the  rigid  array  is  propor- 
tional to  eq.  (40)  , the  average  gain  factor  g for 

. 3/4 

the  free-floating  array  is  proportional  to  K ' . 

This  implies,  for  example,  that  a non-rigid  array  of 

500  elements  has  about  the  same  average  array  gain 

as  a rigid  array  of  100  elements. 

Such  a system  makes  very  inefficient  use  of  its 
data  storage  capability.  It  stores  like  a rigid 
array  system  in  every  memory  element  the  complete 
signal  sequence  received  during  the  long  even  receive 
period.  In  a stationary  array  this  serves  to 
generate  the  coherent  and  focused  odd  insonif ication 
beams  which  carry  not  only  the  spatial  (angular  and 
focus)  information  but  also  the  temporal  (high 
resolution)  range  information.  The  spatial  infor- 
mation is  rendered  useless  by  the  randomization  of 
the  element  positions  and  the  temporal  information 
is  highly  redundant  since  it  appears  in  all  element 
memories  essentially  unchanged.  One  can  argue  that 
the  spatial  information  is  used  to  form  the  inco- 
herently focused  beams  and  are  not  really  rendered 
useless.  However,  the  total  power  at  the  target  is 
only  K times  the  power  of  each  individual  element 
contribution.  One  could,  therefore,  achieve  the  same 
power  level  at  the  target  by  irradiating  it  with  a 
single  wide  angle  projector  capable  of  emitting  K 
times  the  power  of  an  individual  array  element. 
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A possible  system  realization  would  be  to  use 
the  reference  element  as  driver  for  a powerful 
projector  activated  for  the  odd  retransmissions. 

The  reference  element  alone  would  be  equipped  with 
a high-capacity  memory  of  n=2Ate/x  bins  which 

acquires  the  range  information  for  the  targets  in 

\ 

the  field. 

V.  White  Noise 

We  discuss  in  this  section  the  response  of  the 
system  to  additive  white  noise,  that  is,  to  thermal 
noise  at  the  input  of  the  preamplifiers  of  each 
channel.  We  assume  that  the  noise  power  P and  all 
other  pertinent  circuit  parameters  are  identical  in 
all  channels.  The  noise  added  to  the  amplified  and 
filtered  signal  which  is  read  into  the  channel  memory 
during  each  receive  period  is  /P  mn^(t)  where  /P  is 
the  RMS  noise  voltage  and  mn^(t)  a normalized  band- 
limited  noise  process.  The  superscript  m refers  to 
the  mth  receive  period,  the  subscript  k to  the  kth 
channel.  The  processes  ^^(t)  are  statistically 
independent  for  different  m and  k with  a covariance 
Matrix: 

- ■Wskk'a(T>  <57> 

where  a(x)  is  the  autocorrelation  function  of  the 
impulse  response  of  the  filter  following  the  pre- 
amplifier in  each  channel  (see  Fig.  2) . The  <$  are 


the  Kronecker  6: 
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6ae  - < 


1 for  a = 8 , 

0 otherwise 
We  have  assumed  a linear  system.  We  can, 
therefore,  without  loss  of  generality,  consider  the 
response  of  the  system  to  the  noise  input  only  and 
later  add  this  noise  response  to  the  response  of  the 
system  to  a search  pulse  °f  which  we  considered  in 
section  I. 

We  assume  that  the  memory  for  each  element  con- 
tains a sample  °n^(t)  of  the  noise  function.  At 
t — °t  the  transmit  mode  is  enabled  and  the  array 


(58) 


transmits  the  function: 

°fk(t)  = /P  °nk(-t)rect0(t) 
where  rectQ(t)  is  the  rectangular  pulse: 


(59) 


rectQ(t)  = { 


1 for  °ta  < t < 


'3 


(60) 


0 otherwise 

The  target  response  at  the  jth  target  is  according 

to  eq.  (17)  for  point  targets 
o -Si? 


°h • = 


r l D(tjk)  °nk(-t)rect0(t)  = 
3 * 

(7  j /PK 


(61) 


— DUjrnj  (t) 


where 


°nj(t)  = ^ k (D(tjk“tj)°nk(-t).rect0(t) 


(62) 


The  function  “n^ (t)  is  again  a normalized  noise  function. 

The  backscattered  signal  at  the  array  due  to  the 
target  response  °hj  given  by  eq.  (61)  is  in  view  of 
eq.  (19) - (21) 
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(63) 


In  order  to  avoid  unnecessary  complications  we  assume 
that  the  target  ranges  r^  are  sufficiently  different 
so  that  the  backscattered  signals  from  different  tar- 
gets do  not  overlap  in  time.  At  t = °ta  the  receive 
gate  is  activated.  The  signal  ^v^(t)  is  generated  by 
each  element  and  a new  noise  sample  nk(t)  is  added. 
The  combined  received  signal  and  noise  is  amplified, 

O 

filtered  and  read  into  the  memory.  At  t = °tp  = tfl 
the  transmit  mode  is  activated.  The  memory  content 
is  clocked  out  in  reverse  order  and  the  system 
generates  the  first  retransmitted  output  ^fk(t) 
which  is  according  to  eqs.  (22)  , (39a)  and  (63) : 


£k(t) 


ft  D(-°tB)1nk(-t) 


(64) 


The  first  term  on  the  right  side  of  equation  (64) 
represents  the  J retrodirected  beams  generated  from 
the  backscattered  target  excitation  *hj  given  by 
eq.  (61) . The  second  term  is  due  to  the  noise  added 
to  the  received  signal. 

In  order  to  evaluate  the  relative  significance 
of  the  two  contributions  to  the  retransmitted  signal 
*fk,  we  consider  the  resulting  target  excitation  1Pj (t) 
Neglecting  the  side  lobe  contributions  as  irrelevant, 
we  obtain  for  ^Pj (t)  according  to  eqs.  (14)  and  (64): 


r 
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1Pj  (t)  = /PK  g.  D(-tj)°nj  (-t)  ,+  /PKDt-^Vnjlt)  (65) 

where  (t)  is  again  a normalized  noise  function 
defined  analogous  to  "n^  in  equation  (62)  as 


1n  . (t)  =—  1 (Dt ..  -t  ■)  1n.  (-t)  rect.  (t) 

1 /K  k 3k  3 k 1 


with 


rect1(t)  = { 


(66) 


1 for  \ < t < \ 
0 otherwise 


(67) 


The  two  normalized  noise  functions  °nj  and  n^  are 

uncorrelated;  We  have  to  compare,  therefore,  the 

intensities  of  the  two  contributions  to  the  target 

excitation  given  by  equation  (65).  We  find  that 

the  intensity  of  the  first  term  representing  the 

2 

retrodirected  beam  is  (g^)  times  the  intensity  of 
the  second  term.  is  the  round-trip  gain  factor 

defined  in  equation  (39a) . 

We  shall  now  discuss  the  significance  of  this 
result  in  more  detail.  First  we  observe  that  the 
ratio  z between  the  intensities  of  the  backscattered 
and  newly  added  noise  at  the  transmitting  array  is,  as 


follows  from  eq.  (64) : 

z = (g j ) 2/k 


(68) 


whereas  at  the  target  the  corresponding  ratio  is 
according  to  eq.  (65) : 

2j  = (g j ) 2 (69) 

The  ratio  z^/z  = K is  the  array  gain  due  to  the 
coherent  superpositions  at  the  target  of  the  retrans- 
mitted backscattered  noise.  The  system  has  transformed 
this  noise  into  a beamforming  signal  ^s^ft): 
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\<t>  ' 4 ? D(-tj-tjk»°nj(-t> 


(70) 


z is  the  signal  to  noise  ratio  for  this  signal  during 
the  first  receive  period.  If  z is  large  compared 
to  unity  the  target  information  is  already  in  the 
memory  after  the  first  receive  period  and  can  be 
displayed,  for  example,  by  cross-correlating  the 
content  of  one  memory  element  (say  k=0)  with  all 
others  and  displaying  the  correlator  output.  The 
correlator  output  ck(t)  for  the  kth  channel  follows 
from  eqs.  (57),  (62)  and  (64)  as: 

(71) 


ck(t)  - l a{t- (tjk-tj) } 


where  a(t)  is  the  autocorrelation  function  of  the 
filter  impulse  response.  The  function  a(t)  takes 
the  place  of  the  search  pulse  envelope  and  the  dis- 
play is  equivalent  to  that  discussed  above  in  the 
noise-free  case.  If  z is  small  compared  to  unity 


but  gj  large  compared  to  1,  one  obtains  an  improved 


signal  to  noise  ratio  at  higher  iterations.  Using 


eqs.  (65)  and  (19-24)  we  find  for  f. 


k 

f.-  = / xr  l D(-tj-tjk)  (g_j °n j+^n j ) g j + 


2 - = /I 

'k  v K 


(72) 


/P  D(-1tB)1nk(-t) 


We  see  that  the  signal  portion  is  changed  through  the 


incoherent  addition  of  n^ (t) /g^  to  “n^ (t) . The  new 


signal  to  noise  ratio  z is  now 


2z  = g1(l+l/g2)/K 


(73) 


r 
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which  is  improved  by  the  factor, of  g^  compared  to 

the  previous  iteration.  It  is  obvious  from  the 

foregoing  that  every  additional  iteration  contributes 
2 

a factor  g^  to  the  signal  to  noise  ratio. 
nz  = g.2n  Z gr2(V"1)/(K  n£1  g.-2^"1*  ~ (g,2n)/K  (74) 

~ j 


v=l 


v=l 


Equation  (74)  seems  to  imply  that  only  the  round-trip 
gain  g^  of  the  system  determines  the  detectability  of 
a target  independent  of  the  amount  of  the  additive 
noise  in  the  receive  circuitry.  This,  however,  is  not 
true  since  the  total  amplifications  factor  y which 
appears  in  the  expression  for  g^  is  limited  in  any 
realizable  system  by  the  noise  power  Pn»  If  we 
assume,  for  example,  a maximum  rms  voltage  Vmax  at 
the  array  transducers  during  transmission  and  an 
input  impedance  of  at  the  preamplifier  the  maximum 
total  amplification  factor  y is 

* S Vmax//S?V  <75> 

For  a thermal  noise  limited  system  with  Pn  = kTB, 


for  example,  one  obtains 


Y ± Vmax//*1*™  = 1*5  • 10 


v 

10  max 


(76) 


/RjB 


B is  the  bandwidth  of  the  system.  Equation  (76) 
gives,  for  example,  for  Vmax  « 10  KV,  R = 100  fi  and 
B = 10  KHz  a limit  for  the  voltage  amplification 
factor  of  1.5  • lO**  . If  we  introduce  the  maximum 
value  of  y into  the  equation  (39a)  for  we  see  that 
the  maximum  realizable  round-trip  gain  now  depends  on 
the  noise  power  of  the  system: 
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ctflica. 


9jSV 


(RlPn> 


where  a and  3 are  transducer  constants  defined  in 
eqs.  ( 9)  and  (20).  If  gj  approaches  unity,  no 
improvement  in  signal  to  noise  is  obtained  by  higher 
iterations . 

We  shall  now  briefly  consider  the  case  of  a 
search  pulse  together  with  additive  noise  in  the 
system.  The  search  pulse  appears  for  the  first 
time  as  a third  term  in  the  function  1f^(t) 
equation  (64)  . In  order  to  have  the  array  response 
dominated  by  the  search  pulse  rather  than  the  quasi 
random  noise  function  ^sk(t)  we  require  that  the 
pulse  power  P is  very  large  compared  to  the  noise 
power  Pn  and  that  the  round-yrip  gain  g^  is  large 
compared  to  unity.  We  then  obtain  a signal  to  noise 
ratio  for  one  target: 

p g-2 

, * _iT3.. 


which  is  the  equivalent  of  the  case  for  the  quasi 
random  search  function  given  in  equation  (69) 


above . 


VI.  Reverberation  Noise 


In  this  section  the  effect  of  volume  backscatter 
(reverberation)  on  the  system  performance  is  consid- 
ered. Since  the  system  is  linear,  one  can  consider 
again  as  in  the  white  noise  case  the  response  of  the 
system  to  the  noise  only.  It  will  be  found  that  the 
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reverberation  received  in  the  first  receive  period 
is  for  a wide  class  of  systems  practically  uncor- 
related between  different  array  elements  (channels) 
so  that  appreciable  array  gain  can  be  achieved  for 
the  first  retransmit  cycle.  For  higher  iterations, 
however,  the  resulting  reverberation  noise  is  highly 
correlated  between  channels  and  no  further  array  gain 
in  signal  to  noise  ratio  is  possible  with  higher 
iterations. 

The  reverberation  environment 

The  following  analysis  is  based  on  a statisti- 
cally uniform  distribution  of  scatterers  throughout 
the  search  volume.  The  probability  to  find  one 
scatterer  in  a given  volume  element  is  proportional 
to  this  volume  element  and  proportional  to  the 
density  p of  the  scatterers.  One  has  therefore: 

Probability  (one  scatterer  in  dV)  = p*dV  (77) 

Each  scatterer  is  characterized  by  a scatterstrength 
Sj  and  a velocity  v j . The  scatterstrength  Sj  is  a 
complex  random  variable  of  zero  mean  with  the  real 
and  imaginary  parts  jointly  gaussian: 

E{S • } = 0 

3 (78) 

E{S  jS  j i } - 

Vj  is  the  relative  velocity  of  the  jth  scatterer 
relative  to  the  array,  that  is,  the  velocity  compo- 
nent which  determines  the  doppler  shift  of  the  back- 
scattered  signal.  Vj  is  again  a random  variable  of 
zero  mean  which  is  also  assumed  to  be  gaussian: 
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E{ Vj } = 0 

B(viV  - S33' 


6^,  is  the  Kronecker  <$  and  E{x}  is  the  expectation 
of  the  random  variable  x. 


The  reverberation  noise 

We  consider  first  the  signal  received  during 

ci  B 

the  first  receive  period  °t  < t < °tp  caused  by 
the  backscatter  of  the  search  pulse  due  to  the 
distribution  of  point  scatterers  introduced  above. 

One  can  treat  the  scatterers  in  the  same  manner  as 
the  targets  were  treated  in  Section  II  eq.  (19) . 

One  obtains  for  the  pressure  distribution  °pk  at 
the  kth  element: 

nk(t)  * °pk(t)  = 

S. 

E — exp{-iio0(tj+  (tjk+ (tjk-t)  (1+Vj/c)  }°f  (t“tjk-tj)  (80) 

In  eq.  (80)  the  search  pulse  °f(t)  has  been  written  as: 

iu  t 

°f(t)  = e °f(t)  (81) 

where  °f  (t)  is  the  pulse  envelope  and  u)Q  the  center 
frequency  of  the  search  pulse.  The  factor  (l+v^/c) 
in  the  exponent  represents  the  doppler  shift  of  the 
backscattered  signal.  The  doppler  shift  is  signifi- 
cant only  in  the  fast  varying  exponential  but  negli- 
gible in  the  envelope  function. 
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The  noise  covariance 

The  effect  of  the  reverberation  noise  nk(t) 
eq.  (80)  on  the  system  performance  is  governed  by 
the  covariance  function  Rkk,  (t/t)  : 

Rj^.ttjT)  = E{n*,nk(t)  *n*t  (t)  (t+x)  }/E{  |nk(°ta)  |2} 
The  expectatation  E{nk(t) nk, (t+x) } is: 
E{nk(t)nk(t+x) } = 

v . 

E{  E exp{iw0 (t jk~tjk , +t)  }E{exp{iu>0  ^(tj^t^ 


(82) 


,+t)  } 
(83) 


The  double  sum  which  results  if  one  uses  eq.  (80)  to 
form  the  product  nk (t) nk, (t+x)  reduces  by  virtue  of 
eq.  (78)  to  the  single  sum  given  in  eq.  (83).  Since 
Vj  is  gaussian  the  expectation  of  the  doppler  factor 
under  the  sum  in  eq.  (83)  becomes: 

E{exp<iu0  tJk,+T>}  = (84) 

exp{-[w0(v/c) (t-k-tjk,+T) ] 2/B) 

In  order  to  make  the  evaluation  of  the  expec- 
tation of  the  sum  over  all  scatter  contributions 
tractable,  we  consider  only  the  case  of  a line  array 
of  equidistant  elements  and  introduce  a spherical 
coordinate  system  r,  4> , 0 which  has  its  origin  at  the 
reference  element  k = 0 of  the  array,  and  its  polar 
axis  along  the  array  axis.  The  polar  angle  0 is 
defined  such  that  0=0  represents  the  equatorial 
plane  of  the  coordinate  system. 


, 
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The  time  delay  parameters  t.,  in  eq.  (83)  depend 

3 K 

on  the  relative  position  of  the  jth  scatterer  and  the 
kth  array  element: 

V ■ ' ?kl/c  (85) 

where  r^  and  r^  are  the  position  vectors  of  the  jth 
scatterer  and  the  kth  array  element.  In  our  coor- 
dinate system  | r ^ - r^ | becomes  in  first  order  in  (d/r^) 

' . -v  . 

(86) 


'rj  “ rk' 


= Tj  - d k sin  0 


where  d is  the  distance  between  array  elements. 

Introducing  eqs.  (84)  and  (86)  into  eq.  (83) 
one  obtains  for  the  sum  in  eq.  (83) : 

Z = 

j 

1 9 Ire  ^ k * ^ 

o Z — J exp{i0Jo(s+x)-(  (s+T)0Jov/c)  /8}°f  (u.+^fq^-)  °f  (u-j+£TTk  +x) 
j r. 

(87) 

where 


s = (k’-k)d  sin0/c 
2r  j 
c 


Uj  “ * - 


(88) 


To  evaluate  the  expectation  of  the  sum  over  all 
scatter  contributions  eq.  (87)  , one  replaces  the  sum 
over  all  scatterers  by  the  sum  over  all  volume 
elements  which  make  up  the  acceptance  cone  of  the 
system  and  multiplies  each  term  of  the  sum  with  the 
probability  of  finding  a scatterer  in  the  corre- 
sponding volume  element.  One  can  always  make  the 
volume  elements  so  small  that  the  probability  becomes 
negligible  to  find  more  than  one  scatterer  in  the 
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volume  element.  The  occurrence,  of  a scatterer  in 
a volume  element  is  then  a random  variable  which  has 
as  its  range  only  the  values  zero  and  one.  The 
probability  that  the  value  is  one  is  according  to 
eq.  (77)  equal  to  pdV.  The  expectation  of  the  sum, 
therefore,  is  equal  to  the  integral 


E{£}  = 


16pa 

_4 


(t-u) 


-4 


Search 

Volume 


exp{  iu)Q  ( s+x ) - [ wQ  ( s+t  ) ^]  2/8  } 


(89) 


* °5(u  + l^k)of*(u  + jrr=k  + T>dv 

The  coordinate  of  the  jth  scatterer  has  been 
replaced  by  the  coordinate  r of  the  corresponding 
volume  element.  The  new  variable  u replaces  u^  of 
eq.  (88) 

u = t - (90) 

and  the  volume  element  dV  is  given  by: 

2 c4  2 

dV  = r ar  cosO  dQ  d<J>  = § ^ (t-u)  du  ds  d<f>  (91) 

The  integral  is  taken  over  the  acceptance  cone  of  the 
system. 

The  acceptance  cone  is  defined  by  the  radiation 
pattern  of  the  individual  array  elements  and  by  the 
range  R from  which  echoes  are  accepted  during  the  first 
receive  period. 

It  is  sufficient  for  the  purpose  of  the  present 
analysis  to  identify  the  acceptance  cone  with  the  main 
lobe  of  the  radiation  pattern  of  a quadratic  array 


element  of  side  a: 
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+ s "Wx  - arctg  k 


0 < 0 = arctg 

max  2a 


for  a > X/2 


(92) 


The  acceptance  range  is 


(93) 


R . - A - T c < r < R + A + T c 

min  o - - max  o 

where  A = Kd  is  the  array  aperture. 

The  integral  eg.  (89)  can  be  simplified  by 

2 

neglecting  the  doppler  term  exp{- [toQ  (s+t ) v/c]  /8}. 
According  to  eq.  (88)  the  term  w0(s+x)v/c  is  limited  by 
ioQ  (s+t)  v/c  < 2ir  ( j + tv)  v/c  (94) 

Since  v/c  is  of  the  order  of  10  ^ , one  can  neglect 
the  doppler  shift  even  for  very  large  arrays  and 
long  search  pulses.  Introducing  the  volume  element 
eq.  (91)  into  eq.  (89)  and  integrating  over  <J>,  one 
obtains : 

4<(>  o* p 
m 


u 


max 


E{ E j } = 


(k'-k)dt' 


I 1 


'max  iu)0(s+T) 
e 


u 


min 


-|S|  f1  + W 

1 'max 


(95) 


• °f(u  + ^T.ky)  °f  * (U  + + T ) duds 

where  u_  and  u_.  and  Is I are  given  by 

max  mm  1 max ' 


u. 


max 


t Rmin  A 
2 c c ^o 


(96) 


u . = t/2  - J™*  - A/c  - T 

mm  c 


| (k'-k) |d  sin  0 


max 


~ |k'-k|  = I k * — k | vd/2a 


1 max  1 c 

We  shall  now  show  that  one  can  extend  the 
integration  over  u from  -00  to  +00  without  changing 
the  value  of  the  integral  in  eq.  (95) . 


46 


The  time  interval  of  interest  is  the  first 
cc  6 

receive  period  °t  < t < °t  . It  is  given  by: 


A _ t ^ Rmax  , A 
2c  1 c 2 


• + T 

2c  o 


It  follows  from  eqs.  (96)  and  (97)  that  during  the 

first  receive  period  u . is  always  smaller  than 

min 

-A/2c  and  umax  always  larger  than  A/2c  + tq.  The 
— ck 

pulse  envelope  °f(u  + _^) , however,  and,  therefore, 

the  integrand  in  eq.  (95)  is  unequal  to  zero  only  for 

o < u + < Tc  (98) 

Since  by  definition  l^rz^I  is  smaller  than  A/2c,  one 
obtains  as  the  condition  for  a non-vanishing  integrand 
- A/2  < u < A/2  + xo  (99) 

It  follows,  therefore,  that  the  domain  of  integration 
over  u:  u . „ < u < u,„„  always  includes  the  region  where 
the  integrand  is  unequal  zero  for  any  choice  of  t 

within  the  first  receive  period.  The  integral  is, 

therefore,  independent  of  t and  does  not  change  its 
value  if  the  boundaries  are  extended  to  ± °°. 


We  observe  further  that  the  lower  bound  for  t 
(eq.  97)  is  R^.  - j , and  that  the  upper  bound  of  u 

for  which  the  integrand  is  not  equal  to  zero  is: 

u = + To  (100) 

one  obtains,  therefore,  an  upper  bound  for  u/t  by: 

A+2t  c . , 

U/fc  < 2R  - + 0(r^— ) (101) 

2 min  min 

Since  both  A and  tqc  are  small  compared  to  for  all 

practical  systems,  we  can  in  good  approximation  neglect 
the  term  (1  + ^) 2 in  the  integral  of  eq.  (95)  and  obtain: 
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E(V  ■ fkr 


, + s 

+00  1 max  1 


I I 

-00  - 


exp{iu)o  (s+x)  } 


max 


(102) 


• °f  (u  + °f  (u  + + s+x)duds 


Hmop 


where  C = — g — . The  argument  of  °f*  was  obtained 

• V • c Vq 

by  observing  that  + s . We  introduce 

the  new  variables 


sk 


v = u + Rt  _R  and  x = s + x 


(103) 


and  obtain 


E{E  } = 


+00 

C [ | f 1 2dv  Is  I 

J 1 1 1 max 1 

-00  t 


(k'-k)  t" 


- 1 s 


f 1W  X 

A (x) e ° dx  . (104) 


max 


where  A(x)  is  the  normalized  autocorrelation  function 


of  the  pulse  envelope: 

+ 0O 


+oo  +oo 

A (x)  = J °f (v) °f * (v+x)dv/J  |°f(v)|2dv  (105) 

— 00  — oo 

To  evaluate  RkR, {t, (t+x) } as  defined  in  eq.  (82),  we 
have  to  determine  the  function  E{  E } for  t = °t  , x = 0 

j “ 

and  k = k' : 

+°°  2 I s 

E{|nk(«ta)|2  = E{E } = c (°ta)"2  f | f 1 2dv  lim 

3 t=t  ( } 


x=0 

k'=k 


(106) 


+00 


C(°t 


V2  j |f| 


"|2dv  . | v 

a 


Combining  eqs.  (104)  and  (106)  we  finally  obtain  for 
the  covariance  function  RRR, (t,t+x) : 
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(°t  ) .:  ' max' 

Rkk,(t't+T)  - 77^—72 

2(smax)t  -|s  | 

1 max  1 


10)  x 

A(x)e  dx  (107) 


It  is  sufficient  for  the  following  to  suppress 
the  time  dependence  in  .R^^, {t, (t+r) } and  consider 
only  its  maximum  R^,  (°ta,T)  which  we  denote  as  R^i  (t) 
The  integral  in  eq.  (107)  which  is  independent  of  t' 
can  be  formally  integrated.  To  do  so  we  rewrite 
eq.  (96)  as: 


+ 00 


10)oX 


Rj^i  (t)  = — - J rect  (x-x)rect  (x)A(x)e  w dx  (108) 

Is  I J I S I 


max  1 


max  1 


The  rectangular  function  rect^(x)  is  defined  as 


rect  (x)  = { 


1 for  |x|  < y 


(109) 


0 otherwise 


The  first  function  rect 


(x-t)  in  eq.  (107)  was 


max 

introduced  to  change  the  limits  of  the  integral;  the 

second  rectangular  function  rectT  (x)  was  introduced 

o 

to  set  in  evidence  the  fact  that  the  autocorrelation 
function  F(x)  of  the  pulse  envelope  is  zero  for  |x| 
greater  than  the  search  pulse  duration 


A (x)  = rect  (x)  A (x) 


(110) 


a. ) 


To  integrate  eq.  (108)  we  consider  two  cases: 


> T. 


and 


b.)  |s_  | < t 


max'  o ■ max'  o 

The  first  case  corresponds  to  widely  spaced  elements. 
We  shall  show  that  in  this  case  the  correlation  function 


Rkk' is  very  small. 
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To  simplify  the  integral  eg.  (108)  we  consider 
that  the  product  of  two  rectangular  functions  can 
always  be  reduced  to  a single  rectangular  function. 
We  have: 


. rectx  (x) 
for  | t | < 


S < T 

max 1 - o 


rectg  (x-x)rect  (x) 
max  o 


recta(x-8)  (111) 

for  Is  I -T  < T < |s  I + x, 
1 max 1 o 1 max 1 ( 


for  t > | s I > T 
1 1 1 max 1 o 


with 


“ " flsmaxl+To-|Tl)/2 
6 " (To-lsmaJ  + M)/2 

We  need  to  consider  only  the  first  case  which  includes 
the  important  domain  around  x = 0.  Because  of 
eqs.  (110)  and  (112)  the  integral  eq.  (108)  reduces 
to  the  power  spectrum  F(to)  of  the  pulse  envelope: 


with 


*kk'(T> 


2i^~  F(uo)  for  M < 

max 


max  1 


Tc  (112) 


F<%) 


+00 

I 


—00 


ico  X 

A(x)e  dx 


Note  that  (t)  does  not  depend  on  x as  long  as 

lx  I is  smaller  than  |s  I - x 

The  power  spectrum  of  the  pulse  envelope  at  the 
carrier  frequency  w = ojo  is  very  small  even  for  short 
pulses.  We  consider  as  an  example  a pulse  with  a 
gaussian  envelope  and  a time  constant  xQ  = 3/v.  Its 
power  spectrum  F(<*>0)  becomes: 
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„ , . 3/tT 

F(“o>  = 2^r  e 


(113) 


which  is  zero  for  all  practical  purposes. 

We  consider  now  the  second  case  which  corresponds 
to  closely-spaced  elements.  Here  once  again 
correlation  between  elements  is  relatively  small. 

The  product  of  the  two  rectangular  functions  in 
the  integral  eg.  (109)  reduces  in  this  case  to  the 


rectangular  function: 


rectg  (x-x) 
max 


, , , , . for  |t|  < x0  < |smax| 

rect  (x-r)rectT  (x)  = 

smax  To  recta(x-3) 

for  t —Is  I < x < x + 1 s 

o 1 max1  o 1 max 1 


for  t > x + s I 
o 1 max 


with 


“ ■ (IsimxI+to-M}/2 

e - (vlTl-si»ax,/2 


(114) 


We  consider  again  the  first  domain  of  x in  eq.  (114) 

and  find  from  eq.  (109) : 
iw  x 

Rjck1(T)  = e {sinc(smaxu)*F(u)U  for  M < To  " isml  (115) 


The  convolution  of  the  sine  function  with  the  narrower 
pulse-like  power  spectrum  F(gj)  leaves  the  sine  function 
virtually  unchanged  so  that  we  have  in  good  approxi- 
mations for  small  values  of  x: 

Rjck*  (x)  * sinc^osmax)  = sinc(^^U^)  (116) 
For  a dense  array  with  d = a/2  the  correlation  of  the 
noise  between  neighboring  channels  is  zero  in  our  approx- 
imations and  is  not  very  high  for  non-integral  d/a,  which 
is  by  definition  larger  than  unity. 
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Any  coupling,  however,  is  undesirable  since  the 
re-emitted  reverberation  noise  from  correlated 
channels  adds  coherently  at  targets  of  interest  and 
one  can,  therefore,  not  get  the  full  array  gain  in 
signal  to  noise  ratio  during  the  following  iteration. 
One  can,  however,  reduce  the  correlation  even  farther 
by  increasing  the  element  spacing.  A more  effective 
decoupling  of  the  noise  in  neighboring  channels  is 
obtained  by  increasing  both  the  array  spacing  and 
the  search  pulse  bandwidth  so  that  even  for  adjacent 
channels  the  condition  | smax I > tq  discussed  under  a.) 
holds.  This  requires  an  element  spacing: 

d > a(xov)  = n*a  (117) 

where  n is  the  number  of  oscillations  in  the  pulse 
and  a the  element  diameter.  With  such  an  element 
spacing  the  reverberation  noise  in  all  channels  is 
virtually  uncorrelated  and  we  can  treat  the  noise  as 
far  as  its  interaction  with  desired  targets  is 
concerned  in  the  same  manner  as  we  treated  the  white 
noise  in  Section  3. 

We  have  to  consider,  however,  also  the  inter- 
action of  the  retransmitted  noise  with  the  field  of 
undesired  scatterers.  If  the  scatterers  were 
stationary,  the  time  inverted  retransmitted  noise 
signal  would  generate  as  we  have  see®  in  Section  1 
a strong  pulse-like  backscattered  signal  proportional 
to  the  autocorrelation  function  of  the  scatter 


4 


distributions.  The  motion  of  the  scatterers,  however. 
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during  the  time  interval  between  transmissions 
reduces  this  correlation  by  a factor 

,8’vRmax, 


-(2l|Vrto  . lt  ,)2 
e p a e 


-) 


n (118) 
This  factor  had  been  given  in  eq.  (84)  . It  is 
in  the  present  case  zero  for  all  practical  purposes. 

One  can,  therefore,  proceed  with  the  first  iteration 
in  the  same  manner  as  was  done  in  the  white  noise  case 
discussed  above.  With  the  result,  that  one  obtains 
almost  the  full  array  gain  of  K in  signal  to  noise 
ratio.  The  difference  between  the  reverberation 
noise  and  the  white  noise  case  becomes  apparent  at 
the  second  system  cycle  which  begins  with  the 
second  retransmission.  If  targets  are  present,  the 
system  now  transmits  essentially  insonif ication  beams 
focused  on  the  targets.  To  simplify  the  issue  we 
assume  the  presence  of  only  one  target,  therefore,  of 
one  transmitted  beam.  The  reverberation  noise  due 
to  the  retransmitted  beam  (if  we  neglect  the  side 
lobe  contribution  which  was  justified  in  Section  2) 
is  very  similar  to  that  calculated  above  with  the 
important  exception  that  the  acceptance  cone  is  now 
determined  by  the  aperture  of  the  array  A and  not 
the  aperture  of  a single  element  a.  0 therefore, 

111(3.  Jv 


becomes  now 
0 

max 

and  therefore 


. X X 
= arct9  A = Kef 


(119) 
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Is  | = (k~k*  )d0max  .=  (k'-k)  = k'-k 

■ v I n K » \)  V\> 


(120) 


Equation  (123)  shows  that  |sm|  is  even  for  the 
most  distant  elements  much  smaller  than  xQ  . 

We  obtain,  therefore,  in  analogy  to  eq.  (119) 
for  the  covariance  function  of  the  reverberation 
noise  generated  in  the  second  system  cycle 

lRkk'(T)  s = sine  (iJiz^lLL) 

Since  |k-k' | is  smaller  than  K/2  the  noise  in  all 
channels  is  highly  correlated  and  further  iteration 
does  not  increase  the  signal  to  noise  ratio. 

This  behavior  of  the  system  is  demonstrated  in 
the  next  section  by  computer  model  studies. 

VII.  Computer  Simulation  of  System  Operations 
in  Reverberating  Environment 


A computer  model  was  constructed  to  study  the 
behavior  of  the  time  inversion  system  in  a reverber- 
ation environment.  In  this  section  the  model  is 
described,  the  program  listing  given  and  results 
of  the  study  presented.  ^ . 

The  Model 


The  model  calculations  are  based  on  a linear 
receive  transmit  array  which  operates  against  a two 
dimensional  target  field.  A cartesian  coordinate 
system  x,y  is  used  to  describe  the  relative  positions 
of  array  and  targets.  The  array  is  located  on  the  y 
axis  of  the  coordinate  system.  It  consists  of  K = 2n 
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equally  spaced  elements.  The  distance  d between 

elements  is  equal  to  the  wavelength  A of  the  carrier 

' . 1 

frequency  v.  The  unit  of  time  is  the  period  T = — 

of  the  carrier  frequency  and  the  unit  of  length  the 

wavelength  A.  This  convention  reduces  the  frequency 

v,  the  wavelength  A and  the  signal  velocity  c to  unity. 

The  array  elements  are  numbered  from  - ( j - 1)  < k < K/2 

with  the  element  k = 0 coinciding  with  the  origin  of 

the  coordinate  system.  The  target  field  consists  of 

a single  point  Target  of  Strength  S surrounded  by  a 

random  distribution  of  scatterers.  The  target  is 

located  on  the  x-axis  at  x = RQ.  The  scatterers  are 

distributed  over  the  rectangular  area  RQ/2  < x < RQ  ; 

-R  < y < + R , which  is  centered  at  the  target, 
o — — o 

The  total  number  of  scatterers  is  N.  The  position 

> 

vectors  r^  = {x ^ ; y^ } , j = 1,2  . N of  the  scatterers 
are  generated  by  the  computer  as  part  of  the  simulation 
program. 

The  scatter  strength  a ^ of  the  scatterers  which 

appears  always  together  with  the  range  r^  of  the 

2 

scatterer  as  (o^/r^  ) is  replaced  in  the  following  by 

2 2 
its  mean  <aj/rj  > anc*  normalized  to  1/RQ  • The 

normalization  is  possible  since  the  significant  param- 
eter in  the  simulation  is  the  relative  scatter  strength  S 


of  the  actual  target 
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The  Simulation  of  the  System  Operations 


The  target  acquisition  process  begins  as  in  all 

previous  examples  with  the  emission  of  a gaussian 

search  pulse,  which  is  because  of  the  conventions 

adopted  above: 

°f(t)  = exp{i2irt  - t^/x} 

This  pulse  is  scattered  back  by  all  scatterers 

in  the  field  and  by  the  target  proper.  Taking  into 

account  the  normalization  of  the  individual  scatter 

cross  sections,  one  obtains  for  the  first  backscattered 

signal : 

1 N 

x£(t)  ~ { fitt-tj-t jk)  + S • 6(t-t0-t0k)}  * °f(t)  (121) 

with 


V = rjk  " <zj2  + <**-*> 2> 1/2 


t ■ = t . and  t = R 
j jo  o o 

The  function  ( t ) has  to  be  sampled  at  intervals 

At  = where  B is  the  bandwidth  of  the  search  pulse: 


B = 2 
x 

The  sampling  interval  is  therefore: 
At  = x/4 


(122) 


(123) 


The  time  interval  over  which  Xk(t)  is  to  be  observed 

8 a 

is  the  first  receive  period:  °t  - °t  . In  our  model 
because  of  the  convention  c = 1 we  have: 

•t6  - = 2(Rnax  - R^)  + T + K **  “ 4RQ  (124) 

The  required  minimum  number  of  samples  M to  describe 


the  function  X£(t)  is  therefore: 

16  R 

M At T~ 


(125) 
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It  is  more  efficient  to  perform  the  actual 
calculation  in  the  Fourier  domain.  The  Fourier 
transform  X^(v)  of  x£(t)  is  given  by: 


x£(v)  * { 2 e 

j=l 


N -i2ir  v (r  .+r  ..) 


j + 


-i27rv(r0+r0k) 


- (tt  (v-1)  t ) 


The  sampling  intervals  Av  are  given  by: 


(126) 


Av  = 


2(°t0_ota) 


(127) 


and  since  M samples  have  to  be  taken,  v ranges  between 


1 - M/2  Av  < v < 1 + M/2  Av 
1--<v<1  + - 

T T 


(128) 


The  frequency  band  is  centered  around  the  carrier 
frequency  which  according  to  the  adopted  convention 
is  equal  to  1. 

The  calculation  of  implies  the  calculations 
of  K x M sums  over  the  echoes  from  the  N scatterers. 
The  actual  results  shown  later  are  based  on  1000 
randomly  distributed  scatterers. 

In  order  to  reduce  the  amount  of  calculation. 


we  replace  rj~  by  a linear  approximation 


rjk  = 


2 — 2 

+ (yj-kr  - r.  - 


(129) 


and  replace  the  sum  over  j in  eq.  (126)  which  has 
in  our  model  1000  terms  by  a sum  over  100  equally 

spaced  phase  terms  in  the  following  manner: 

y • 100  i2ir 

N i2ir{v(2r.  + k -^) } = Z a(v,k)0  • e u 

Z e 3 rj  JL=1  * 


(130) 


I 


where  a(v,k)^  is  the  number  of  occurrences  that  the 

yi 

fractional  part  of  v • (2r^  + k p*4-)  falls  between 

£,-1  1 3 

and  £./100.  The  Xk(v)  calculated  in  our  program 

is  therefore: 


i 100  i2rrvR^ 

x£(v)  = { E a(v,k)p  el2lT(jl"*5)/100  + Se  °} 
£=1 


-{  (V-l)  *TTt}‘ 


(131) 


The  computation  of  x£(v)  is  executed  in  two 
steps.  The  first  step  evaluates  X^(v)  witn  no  target 
present  that  is  for  S = 0.  The  result  is  then  put  on 
tape  (identified  as  tape  16  in  the  program) . The 
second  step  is  to  combine  this  with  a target  response. 
This  partition  has  been  made  to  facilitate  experimenting 
with  different  target  sizes,  locations  and  numbers. 

The  program  listings  are  given  in  Tables  1 and  2. 
The  inputs  to  the  program  are  the  pulse  duration  t 
identified  as  TW,  the  target  position  RQ,  the  number  of 
array  elements  K,  the  target  strength  S and  the  total 
number  of  scatterers  N. 

The  Second  and  Higher  Iterations 


The  second  iteration  X^(t)  represents  the  return 
from  the  retransmitted  time  inverted  signal  Xj^(t)  an(* 
is  given  by: 

x£(t)  - tu  + s/6(t-tok-tok,)> 

k ] 

• *x]^|')  (-t) 


(132) 


The  term  Yj  in  the  time  delay  operator  of  eq.  (127) 

represents  an  approximation  to  the  effect  of  the 

motion  of  the  scatterers.  Yj  is  a random  number 

which  represents  the  displacement  of  the  scatterer 

which  occurred  between  successive  pings.  Yj  is 

gaussian  with  a variance  of  2u  generated  in  the 

program.  Introducing  the  same  approximations  and 

conventions  as  in  evaluating  we  obtain  for  the 

2 

Fourier  transform  of  X^(v): 


X(-v)  * E {E  e 
k k'  j 


-i2Trv(2r‘  + y. (k-k')/r,  + y.) 


3 J3 


3 '3‘ 


-i47rvr. 


+ Se 


(133) 


• x£, (v) 

Here  again  we  approximate  the  sum  over  j by: 


-2tt{v  (2r  • + y.(k-k')/r  + y.} 

E e 33  3 = E(vk)  = 

3 


100 

E a(vK)e 
1=1 


— 2ir  (x,-.  5) /100 


where  k = (k-k1).  k can  assume  2 K different  values 

2 

for  which  E(vk)  is  evaluated.  Xv(-v)  then  becomes: 

, k+K/2  -2ttvR  . 

X^(v)  - E (E(vk)  + Se  )Xk-K(v) 

<=k-| 

2 

The  calculation  of  Xk(-v)  is  again  executed  in  two 

steps  - first  for  S = 0: 

, k+K/2  , 

{Xj"(-v)  } = E E(vK)x£  (v) 

s=°  r-v  K 
K_k2 
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and  then  for  the  return  with  target 

2 - k+K/2  -4irvR  , 

X(-v)  = {XjM-v)}  + E Se  • x£  (v) 

k S=0  K=k-K/2 

The  program  listings  are  given  in  Tables  3 and  4. 

The  Fourier  transform  of  the  first  and  second  and  high 

1 2 

iteration  returns  Xk(v)  and  X^(v)  were  transformed 

back  into  the  time  domain  by  the  Fourier  transform 

and  plotting  routine  given  in  Table  5. 

To  compute  the  third  iteration  the  return- 
2 

functions  Xk(v),  which  is  given  by  the  program  in 

the  frequency  domain,  has  to  be  transformed  back 

2 

into  the  time  domain  X^(t)  and  multiplied  with  a 

window  function  to  account  for  the  narrow  odd  receive 

periods  (see  Fig.  1) . This  truncated  set  of  func- 

tions  X^t  is  then  transformed  back  into  the  frequency- 
~~2 

domains  to  give  X^(v). 

The  third  iteration  then  follows  in  complete 

analogy  to  the  second  one  by 

k+K/2  -4ttvR  ~ 

Xj*(v)  * Z (I(-v,K)  + S»e  }X‘(v) 

K <=k-K/2  K 

Results  of  the  Simulations 

The  program  described  above  is  capable  of 
handling  multiple  targets  and  a wide  range  of  search 
pulse  bandwidths. 

Due  to  budgetary  constraints,  it  was  only  possible 


to  consider  a small  system  of  32  elements  operating 
against  a single  target  with  a narrow  band  search  pulse 
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A target  range  of  Rq  = 1000  and  a pulse  width 

of  t = 160  was  chosen.  This  made  it  possible  to 

handle  the  return  with  128  frequency  or  time  samples. 

In  Fig.  6a  the  return  signal  X*(t)  without  a 

target  in  the  field  is  shown  for  all  32  channels. 

2 

In  Fig.  6b  the  first  iteration  return  X^(t)  is 
shown  for  the  case  that  all  scatterers  are  stationary. 

The  return  signal  to  the  reference  element  is  essen- 
tially the  autocorrelation  of  the  reverberation  noise. 

The  signal  in  the  channels  k / 0 gives  a measure  of 
the  noise  correlations  in  separate  channels.  It  can 
be  seen  that  the  correlation  is  still  substantial  in 
the  nearest  neighbor  channels  but  has  practically 
disappeared  in  all  other  channels. 

In  Fig.  6c  the  first  iteration  return  is  shown 

with  the  scatterer  motions  taken  into  account.  The 
2 

second  return  X^Ct)  is  now  essentially  noise-like 
similar  to  Xj^(t)  as  is  to  be  expected. 

In  Fig.  7a  the  return  signal  X^(t)  is  shown  with  a 
target  of  strength  S = 15  present  in  the  field.  The  signal 
still  looks  essentially  noise-like.  The  target  can- 
not be  recognized.  Figure  7b  shows  the  first 
2 

iteration  Xk(t)  with  the  target  in  place.  The  target 
is  now  clearly  visible.  This  is  due  to  the  residual 
array  gain  which  is  smaller  than  the  array  gain  (=K) 
against  uncorrelated  noise,  but  which,  due  to  the  low 
correlation  of  the  noise  in  channels  more  than  two 
elements  apart  (see  Fig.  6b) , is  still  appreciable. 
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VIII.  Conclusion 

A new  concept  of  a retrodirective  self-focussing 
sonar  or  radar  system  has  been  introduced.  The  new 
concept  is  based  on  time  reversal  effected  by  storing 
and  re-emitting  the  received  signals.  The  memory 
requirements  for  its  implementation  are  discussed 
and  it  was  shown  that  state  of  the  art  serial  access 

memories  are  compatible  with  the  system  requirements. 

. 

The  motion  sensitivity  of  the  system  was  analyzed. 
It  was  found  that  the  system  is  compatible  with  an 
aggregate  of  free-floating  sonobouys0  It  was  further 
shown  that  the  motion  effect  on  a system  mounted 
rigidly  on  a moving  platform  can  be  largely  compensated 
by  an  appropriate  shift  in  the  clock-rate  which  governs 
the  data  flow  through  the  memory. 

A computer  model  to  simulate  the  response  of  the 
system  to  reverberation  noise  has  been  developed  and 
completed.  However,  funds  were  insufficient  to  use  the 
model  for  more  than  the  simulation  of  a very  simple 
low  data-rate  case. 


— *J 
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Appendix  I 


Listing  of  the  Program 


to  Calculate  Iterative  Returns  Without  Noise 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
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1. 

2. 

v 3. 

4. 

5. 

6. 

7. 

a. 

9. 
10. 
11. 
12. 
. 13. 
14. 
l*i. 
16. 

17. 

18. 

19. 

20. 
21. 
22. 

23. 

24. 

25. 

26. 
27. 
23. 

29. 

30. 

31. 

32. 
32. 

34. 

35. 

36. 

37. 

38. 

39. 

40. 

41. 

42. 

43. 

44. 

45. 

46. 

47. 
44* 

49. 

50. 

51. 
62... 
53. 
.54. 

--  55. 

56. 

57. 

58. 

59. 

60. 
61. 
62. 

63. 

64. 


PROGRAM  bONAR  ( INPUT . OUTPUT » TAPF.b=INPI  IT  * TAPr 6=0UTP()T ) 
CO-vOf!  X < 5 I »SIGva(s>  # R 10 ( 6 ) » RIK  ( b » SJ  3)  . pJJ ( 5#  5 > 
COMMON  C-JJ  ( 6 » 6 ) > LJJK  ( 5 » b . 32  J » KUJK  ( S . b # 3?  > 

COMMON  Y H ( b 1 2 ) . YH(512I . THAM204H)  »FKC(2(148) 

PEAL  NU.LJJK.KJJK.NU1 
COVVON/COVKFT /TO  ( 1024  1 .NTAIUE 

, COMPLEX  FKC#  TO.  SUM*  TERM. T HAS. SUM? .CONST 7 

!PI=3.141b 

PI2=2#3. 141b  : 

WHITE  l fa#. 3)  i ':'1  7_  \y 

3 FORMA!  (6X»  *1  ARGET  NO  .*»  5X * *SIGMA*  * 10X > *X*  , 1 1 X#  *?.*) 

REAP  (5.11)  NTARG 

11  FORMAT  (12)  

PCI  1 )=1#NTAPG 

FLAP  (5.12)  S10‘.'A(  I)  »X(  I ) >7 

12  rORPA  T ( 3<F16.  U»  5X)  ) . 

WR I TE  ( 6 » 13  > I .SIGMA (I)#X(I)#Z 

13  FOPMAl  (7X,nx.  I2.7X#F_15.6»3X»E15.b»3X»E15.5> 

DO  2 J-1'613 

PIK(T»J)=S0RT(7*Z4(X(I)-XK)**2)  1 

2 XK=XK*1.  . 

PIO( J )=RlKl 1 #257) 

1 CONTINUE 

READ  l b# 12)  T All#  GAMMA  _ 

WRITE (6.141  TAN. GAMMA 

14  FORMA  TAII  = *»F15.7»*  GAMMA  = ♦.Fib. 7) 

. PC  5 1 = 1 .NT APS  - • ...  .. 

CONST=XU)/K10(I> 
nO  6 J=1»N1ARG 

. _ RJJ( t » J)=C0M5r-(XlJ)/RI0( J) 1 - 

6 CONTINUE 
5 CONTINUE 

DO  7 I = 1»NTARG  . — - .. 

DO  P J=1 #NT ARG 

CUJ(I» J)=2*(RI0< I)-RIO(J) > 

fl  CONTINUE  . _ 

7 CONTINUE 
XMAX=0. 

DO  JO  J=1#NTAR6  ....  ._  ...  . _ 

10  XMAX  = XMAX-M  ( SIGMA (J)*GAMMA>/ (RIO (J)*RIO(J) ) )**2 

XMAX=bl 3. 4XMAX 

jf"AyrXMAX«  1,3  ... 

XMAXl  = XMAX/ll)OnO. 

TKJMTN=0 

DO  20  1 = 1#  NT  ARG  

TF(X(I).LE.O.)  21 #22 

21  TK JTFMPn-R IK (1.513)+RIO(I) 

..  GO  TO  23  . - 

22  TKJTri'fJ=-RIK(  1, 1)41(10(1) 

23  IF (TKJTEMp.LT .TKJMIN)  TKJMIN^TK JTEMP 

. 20  CONTINUE  ...  

TKJUAX=P. 
no  30  I=1#N1ARG 

IF ( X ( I ) . GT . 0 . ) 31 » 32 

31  TF ( y ( I ) .GT.256.)  33.34 

33  ITENP=613 

GO  TO  35  1 . 

34  ITFVP=X(I)4257. 

GO  TO  35 

32  !F(X(I) .LT.-256.)  36#37 

36  ITff'psi 

60  TO  35 

.37  TTFMPSX  1 1)4256. 
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65.  . 35  TKJTrMP=-RlK(I*ITEMP>+RIO(I> 

66.  IF(TKJTFMP.GT.TKJMAX)  7KJMAX=TKJTEMP 

67.  30  CONTINUE  . ....  . — 

68.  ■ nELTANU=PI/(TKJ*'AX-TKJMIN+R!2*TAU) 

69.  TMAX=(TKJV«X-TKJMlN4PI2*TAU)/2. 

70.  MT=9*PT/ (TAIHOELTANU) 

C-  CALCULATE  1.06  HASr  2 OF  NUMBER  OF  DATA.  POINTS 

71.  DO  110  UbTARF=l»2n 

72.  IF  (2**NS) AGE.GF .NIT)  GO  TO  120  

73.  110  CONTlNUt 

79.  120  NIT=2**MSTAGE 

75.  MU5TRT=-riIl*PEl  TANU/2.  ... 

C-  I5TEP  TS  THE  INCRF^ENT  FIETWEEN  ARRAY  ELEMENTS 

C-  FOR  V.T'TCH  DATA  IS  TO  PE  CALCULATED 

C-  1STEP  SHOULD  Ft  A POWER  OF  2 * 1 .LE.  ISTEP  .LF.  256  - 

76.  ISTEP=3? 

77.  f.'K=512/ISTEP 

C-  CALCULATE  LOJK.KJJK  - • • 

78.  DO  ?‘>0  J=1 * NTAPG 

79.  no  260  K=1*NK+1 

C-  IK.  IS  THE  ACTUAL  ARRAT  ELEMENT  WHEN  THEY  ARE  NUMBERED  (1*513) 

00.  IK=<K“l)*lSTff>41 

01.  TE'  P1  = TKJMIN4R1K( J* IK)-RIO(J) 

62.  T£MP2=TKJNAX+RTK(  J.  IK)-RIO(  J) 

83.  DO  250  Jl  = l *NTARG 

C-  DO  rOT  USE  IHIS  OH  LEADING  TERMS 

89.  . . JF(J.EO.Jl)  CO  TO  250  

05.  KJJK ( J. Jl *K )= (TEMPI +CJJ(J*Jl > )/ABS(BJJ< J* Jl) ) 

86.  I.JJK ( J* Jl *K)  = (TFf'P2  + CJJ( J* Jl) )/ABS(HJJ(J* Jl ) ) 

C-  KJJK  CAN  BE  NO  LESS  THAU  -256  . ._ 

87.  TF(KJJK(J*Jl*K>.LT.-?56IKJJK(J*Jl*K>=-?66 
C-  LJJK  CAN  PE  NO  MOPE  Than  256 

88.  IF (LJJKt J* J1*K> .CT.256)  L J JK ( J * Jl * K ) =256  

69.  250  CONTlNUt 

90.  200  CONTlNUt 

91.  290  CONTlNUt 

C CORRECUON  FACTOR  SORT  (PI ) +TAU*OELT AMU/  (2*PI > 

92.  pISCR  T=T  AU*UELT  ANU+ .2820998 

93.  ...  TAUS09- 1 I AU*TA())/9.  ■! 

99.  NK=NK/2 

95.  DO  999  K=1*NK+1 

C-  ZERO  CUT  FXC  ._  . _ _ 

96.  DO  950  T1=1*204H 

97.  950  FKC ( 1 1 ) - ( 0 • * 0 , ) 

98.  IK=(K-1  )*ISTEP+1  -.  . 

99.  NU=fllSTRl 

C-  IK1  IS  ACTUAL  ARRAY  ELEMENT  WHEN  NUMBERED  (-256*256) 

100.  IKlsIK-257  • 

C-  CALCULATION  OF  FK(NU) 

C 

._  101..  . no  990  I = 1*NIT  L __ 

102,  SUvr(0.*0.) 

103.  SUM2=(0.*0.) 

109,  TERr=(0.*0.)  . 

105.  no  980  J=1.NTARG 

106,  C0MST=PI2+NU 

...  107.  C0).ST1=P)S0RT*EXP(-NU’*NU*TAUS09) 

103.  A J-  (SIG’'A(J)*GAMMA)/lRJO(J)*HIO(d)  ) 

109.  ro  970  Jl  = l *MTARG 

- 110.  IF(J.EO.vU)  GO  TO  970  

111 . AJl=(SIGMA(Jl H GAMMA )/(R 10 ( Jl ) *KIO ( Jl ) ) 

112.  BJ=ABS(RJJ( J* Jl ) ) 

C-  IF  LJJK  .LT.  KJJK  HASH  TERM  IS  ZERO..  
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113. 

114. 
1155. 
116. 
117. 
llfl. 
110. 
120. 

I 

121. 

122. 

123. 

I 

124. * 

125. 

126. 
127. 
120. 

129. 

130. 

131. 

132. 

133. 

134. 

135. 

136. 

137. 

138. 

139. 
146. 

141. 

142. 

143. 

144. 

145. 

146. 

147. 

148. 

149. 

150. 

151. 
162. 

153. 

154. 
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IFILJJK ( J# J1 #K! .1  T.KJJK ( J» J1 »KI » GO  TO  470 
CONST 3-  ( C'ONST*P  J ) /2 • 

C0tiST4  = LJJMJ».Jl  «K)-K.UK(J*  JliKIM  . ..  ... 

C0f.ST5-SJN(  CONS  T 4 ♦CONST  3 ) / SIN  ( C0NST3  ) 

C0NST6=-tvJJ(  Ji  J1  I ♦ (LJJK  < J»  JI  i K ) + K JJK  ( J • J1  * K 1)  J/2 . 

TEH”::  TERM*  A JI  *CnfiSTS*CFXP<CMPLX  tU.  t -CONST  *C0NST6)  ) 

470  COtITINIIt 

TKjrplO(J) “R I K ( J » IK  ) 

C-  ADO  HArH  TERM  AMO  LF  At)  I MG  TERM  - 

CONST  /=AJ*CFXP(  CHPLX  1 0 . i “CONST  * IK  J)  ) 

SU''?rSU**2*  TERv*CONST7 

SUf-SUW+AJ+613. *C0NST7  .._  

C-  MOW  MULTIPLY  HY  E.XP ( ( -NU*T AU/2 1 **2 1 
"4110  CONTIMUt 

FKCtI  ) = (6UM-*SUM2)  ♦CONST1  ..  ..  

TH/.S(n=6HM2*C0NSTl 

NU“MU+OELTANU 

490  CCNTINUF.  

WHITE <6. bid)  Ik  1 
C-  COMPUTE  FET 

NTAPLE=0  

CALL  FFT IfJSTAGE # 1 • 0 »FKC 1 
CALL  FFT(M5TAGF,1.0.THAS) 

C-  SHIFT  IM  TINE  

C.AIL  5H  ( Nil  »FKC  ) 

CALL  SII(MITiTMAS) 

C-  PRINT  OUT  1 

510  FORMAT ( ///*  ARRAY  ELEMENT  =*»X5/1 
• WRITE (6  » 520) 

520  FP|-MATI/1?X»*T*»10X.TFK(T)»*10X#*HASH*/)  

521  FORMAT (/12X.*Mlj*,inx»*HEAL(FK)*»10X**IMlFK)*/) 
PELTAT=2.*TMAX/NIT 

T=TPAX  . ...  

no  560  t = i * m i r 
YH( I >=CAHS«  THASI 1 ) ) 

. YPI 1 1-CAPS (FKC  (111  . 

V.'HITE  (6*640)  T » YP  ( I ) » YH  ( 1 1 

540  FORMAT ( 1X»F15.6*?E 18.51 

541  FORMAT ( 1X«F15.5»2E15. 51  ..  

T=T-nELTAT 

550  CONTINUE 

CALL  PLOI  (YP»MIT»XMAX)  , I 

CALI.  PLOI  (YH.NIT  * XM AX  1 ) 

CAI  L PLOI <YH»NTT#-1.) 

999  CONTIMUt  1.  .1..- 

STOP 

END  . 
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1.  SUpPOUTTME  EFT  (NST AOE#PSIGH*X ) 

C-  COPYRICHT  1973  PY  J.  PRESTCOTT  AND  R.  J.  JENKINS 

C DEPARTMENT  OF  ELECTRICAL  ENGINEERING  SCIENCE  ' 

C UNIVERSITY  OF  ESSEX  COLCHESTER  ESSEX 

C THIS  SUtIROUTINF  MAY  HOT  RF  USER  FOR  COMMERCIAL  PURPOSES 

C WITHOUT  THE  EXPRESS  WH II  TEN  PERMISSION  OF  THE  ABOVE. 

C — - - - 

C THIS  SUBROUTINE  IMPLEMENTS  a fast  FOURItR  TRANSFORM 

c algorithm  given  or  m.  l.  uh rich  i/eee  thans.  audio  and 

_C  ELFCTRP ACOUSTICS*  JUNE  1909*  170-17?)  I AMD  W.  T.  COCHRAN 

C ET  AL  <HROC.  IEEE*  OCT.  lRf>7*  SSI  lbG4-lb7‘)) . 

C 

c .....  

C THIS  SUBROUTINE  IS  CALLED  BY 

C 

C CALL  FFTUISTAGE.RSIGN)  .- - 

C 

C where  NSTAGE  IS  THE  LOG  (RASE  2)  OF  THE  MUMPER  OF  POINTS  IN 

C The  TRANSFORM#  SO  THAT  FOR  12H  POINTS  NSTAGF=7.  _. 

C The  SIGN  GOVERNS  THE  DIRECTION  OF  THE  TRANSFORM#  AND  IS 

C 

_C  -1.0  FOR  FORWARD  TRANSFORMS  

C +1.0  FOR  INVERSE  TRANSFORMS 

C FIRST  N LOCATIONS  OF  ARRAY  X.  THE  NORMAL  I ZED  FORWARD 

.c  v.  transform  or  the  inverse  transform  will  appear  in  the 
c the  n complex  points  of  the  input  sample  shoulo  hf  in  the 
c 

. c correct  order  in  the  same  n locations.  - 

c 

c thf.  sample  array  x in  comfft  must  have  at  least  2*n  elements 

C and  the  COEFFICIENT  ARRAY  I must  have  at  least  n elements. 

C 

c the  suproutinf  stanos  alone  and  can  he  usEn  for  ant 

C COMnlNAIION  OF  FORWARD  AND  INVERSF  TRANSFORMS  of  any  size • 

C without  outside  attention. 

• c 

2.  COMPLEX  X»W*TX1»TX2»T 

3.  INI  EGER  SRC1#SRC2#DST1»DST2 

4,  CCJ'MON/COMPFT/T  ( 1 02  <4 ) * NT  ADLE 

5.  DIMENSION  X ( 204  4 ) 

< 6,  DATA  NTAHLE/O/ 

_7.  TPI=3.0tATAN(1.0)  

c 

C COMPUTE  COEFFICIENT  TABLE  AND  SET  NTAPLE  TO  PREVENT 

C RECALCULATION  ON  SUBSEQUENT  CALLS  TO  THE  SUBROUTINE. 

C 

8.  IF (NT ABLE. GT .0)  GO  TO  30 

9.  PP=0.0  . 

10.  DO  20  1=1. S12 

11.  T ( n=CMPLX(C.CS(PP)  *SIN(PP)  ) 

_.  12.  — i T(512  + I)=CONJGtT(D)  — 

13.  20  PP=PP-TPI / 1024.0 

14.  NTAPLE=1V24 

15.  30  CONTINUt  ...  ..  — - * - - - 

c 

C SF.T  UP  I HE  LOOP  CONTROL  VARIABLES 
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16. 

17. 

C 

C 

C 


10. 

19. 

20. 


C 

C 

C 

c 
c 

23. 

c 

c 

c 
c 

24. 

25. 

C . 

C 

C 

26. 

27. 

29. 

29. 

C 

C 

C 

30. 

31. 

32. 

C . 

C 

C 

33.  \ 

34. 

35. 

. 36. 

• C 
C 

. 37. 

38. 

39. 

40.  .. 

41. 

42. 

. 43.  __ 

44. 

45.  1 

46.  2 

47.  3 
C 

- C--  - 

C 

48. 

..  _ 49.  ' 

50. 

51.  4 

52. 

53. 


INCRrNTAPLE 

1U0XIT=1 
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FOR  INVtRSF  TRANSFORMS  SELECT  THE  COMPLEX  CONJUGATE  TABLE 


iF(PSlGfJ.GT.O.O)  INnXIT=(MTABLE/2)+l 

N=2*«NSTAGE 

N2=N/2 

JLIMIT  = 1 . _ - 

KL1MIT=N 


SET  DST2  TO  N 10)  FOR  NSTAGE  ODD  1EVEM1--70  GOVERN  THE 

CORRECT  ALTERNATION  OF  COMPUTATION  BETWEEN  HALVES  OF  THE 
ARRAY 

DST?=N.AND. lB192+204ft+S12+128+32+8) 

IF  AN  INVERSE  TRANSFORM . THE  NORMALISATION  LOOP  WILL  BF  

OMITTEO*  SO  CHANGE  THE  ALTERNA I ION  SEOUENCE. 

IF(RSIGN.GT.O.O)  DST2-0ST2. XOR .N  . 

SRC2=0 

OUTER  LOOP  

00  3 1=1 » NSTAGE 

KLlMlT=KLIMIT/2  — 

INCR=INCK/2 
INDX=INnxiT 

ALTFHNAIE  THE  COMPUTATION  RETWEEN  SUBARRAYS  AS  REQUIRFO. 

DSTl=OST2.ANO.H  _ - ...  - - — - - - - 

SRC?=(SRC2.AND, .N0T.DST1) .ANO.N 
PST?=DST1+N2 

MIPPLE  loop  selects  the  coefficient 


CO  2 J=1»JLIMIT 
SRC  lrSRCU 
SRC2-SRC24KL1MTT 
S'=T  (INFIX) 

1 Nf;FR  LUOP 


no  1 K=1»KLIMIT 

SRC1=SRC1+1 

SRC?rSRC241 

DSTlsOSTl+l  

rsT?=usT2+i 

TX1=Y ISRC1 ) 

TX2=W»X(SRC2>  , 

XIPST1 )=IX1+TX2 
X(rST2)=TXl-TX2 

lNnX=If:l.’X'*  JNCP  

JL IMIT=JL1MI t+jlimit 

IF  FORVAHO  TRANSFORM  THEN -NORMALISE- 

IF(RSIGN.GT.O.O)  RETURN 

PP=1 . 0/FLCAT  (N)  

00  4 1 = 1 *N 
X(I)=X(N4I)*PP 
RETURN 
ENO 
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C-  IF 


SUPROUT1NF  PLOT(X.NN.XMAX) 

C-  PLOTS  N POINTS  IN  X ON  PRINTER.  XMAX  IS  X VALl'P  WHICH 
C-  IS  TO  RE  MAXIMUM  ON  PLOT.  IF  XMAX  .LE.  0.  PROGRAM  WILL 
C USE  LARGEST  VALUE  TN  X AS  XMAS 

PIMEMSION  XII ) .CHARI  100) 

INTEGER  ULANK  * PL  * CHAR  . - 

DATA  NCOL/1UO/ 

N=NH 

PL=1H*  ...  — ... 

PLANK=1H 
WRITE  1 6* 1 ) 

C-  IF  XV AX  .LE*  0 MAXIMIZE  X-ARRAT..-IO.  GET.  XMAX 

IF  ( yuAX . (3T  . 0 ) GO  TO  10 
DO  •?.  1 = 1*  N 

IKX(I).OT.XMAX)  XMAX=XII)  

2 CONTINUE 
C-  NOW  FILL  Char  RUFFER 

10  DO  SO  1 = 1 *N  . 

I X= ( X ( I ) /XV AX ) *MCOL 
IF( IX.GT. 11W)  7X=100 

IFITX.LT. 0)  ixso  

DO  ?0  1I=1*MC0L 
?0  CHARI  I 1 >=MLAMK 

IFIIX.ro. 0)  GO  TO  <40 * 

DO  .TO  1 1 = 1 * I X 
CHARI  II) =PL 

30  CCNUNUE  . . 

C-.  PRINT  LINE 

40  IE(IX.tQ.O)  CHAR  1 1 ) = 1HI 

1 FORMA  I (//IX. * — — 

WRITE  16* 100)  CHAR 
50  CONTINUE 

100  FORMAT  (1X.10UA1)  

WRITFI6.1) 

RETURN 
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to  Simulate  Reverberation  Environment 


» 


t 
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TABLE  1 
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] • 

prcoo*  rv  <*  i -to  ( r ,,!-t  1 7 » i i8 » T.AHt.  | y » outpj  r « 

2. 

COfli'iuj.MK  , 2n)».i\(K  » iOUCi)  *u  ( 1000) 

3. 

Ulfc£N$riN  X ( 1000)  »«  (2000)  *F  { 2n  ) , T ( 2n  ) . 

4. 

ulf.cii'.S  T •'  'i  Y ( 2n  ) 

5. 

J >1  &.*,’(  lOO)  * AL  (100) 

6. 

C 0M2 L t * Y « 2 » T Em^ » TK » LX 

7. 

1 f»  f E6££  S , T . 

8. 

U#\ T A r . ►< fi  * f\  » s » 0/  open 

9. 

HIs4.o*MAN(1.0J 

10. 

Pa2*t>J 

11. 

LL  = «*-\ 

12. 

UK=1.U/1  1 70^.0 

13. 

SEE'.  ( 1 > =0.0 

14. 

SEED ( 2 ) = 1 .0 

15. 

CALL  «<A’  3F  {-2000  t3E£0) 

16. 

CmLL  i3F  (2000  .*<*SEED) 

17. 

Du  1C  T = 1 . J 

18. 

H < I ) = -ffl  ♦ 1 0 0 'J®1-  ( 1 > -5  v ./ 

19. 

A ( 1 ) = 2 r.  i ) ■»  S(  ( ! i ♦ £ ) - 1 0 0 w 

20. 

Ao=-<  ( f ) »#:>♦/.  < I ) «<»2 

21. 

K ( 1 ) = I ( Ai , ) 

22. 

u ( I ) = .S  ( T ) / H ( I ) 

23. 

uO  7 .)  = ’••< 

24. 

7 

HK  ( J.  1 1 =w  ( I ) - ( J-lf>)  * j ( I ) 

25. 

IT 

coo  r i j il* 

26. 

DO  1 t = I . 1 0 J 

27. 

1 

LA  ( n rr--^(CHML.\  U»J»H»{.U1*I-.0QS)  ) ) 

c8. 

uu  3 i = i *lL 

29. 

Y ( i ) = l . / I -LL/2) #9r 

30* 

3 

T ( I ) = r < - < - ( l ' 1 «( 1 - F < 1 ) ) * T .• ) * * 2 ) 

31. 

DO  11  <«  = l,r. 

32. 

DO  10  T =1,lL 

33. 

DO  ± . 1 00 

34. 

') 

«L  ( J ) =,i . O 

35. 

DO  5 J* 1 * 4 

36. 

L*=f-  (1  ) * (Oi<  (r-K,j)  >iv  ( j ) ) 

yr. 

r-'r  (0)  ) *10u*l 

38. 

5 

AL  ( 1)  s.M.  ( 1)  ♦ 1 

39. 

Y(1)=(O.O.U.O) 

40. 

Du  o J= 1 » 1 0 0 

41. 

6 

Y ( I ) » Y ( T ) ♦ AL ( J)  #Ea  ( J) 

42. 

Y(i)=YI» >*T<1) 

43. 

Z(f\i\* I ) = Y < I ) 

44. 

l() 

Con  i i N't*- 

45« 

11 

COnT  i-4'ic 

46. 

WDI It  (l  n ( (/  (J.i) »J=i .32) *1=1.128) 

47. 

S 1 OF 

48. 

Enu 
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TABLE  2 

THIS  PAGES  IS  BEST  QUAllTT TBACttCAK* 

PROMCX>PIPUraiISHBI)I»DDC.  


1 . 

PROGRAM  lM6(TAPE17,TAPE:i0*TAPF19.0llTPOT.TAPE6=0UTPUT) 

2. 

COMMOM  Z ( K .»  2n  ) » T ( 2n  ).F(2n  ) 

3. 

D If«'ENf, I ON  Y(2P  ) 

4. 

C09PL rX  Z.Y 

5. 

INltGF.P  S*TW 

6 • 

DATA  TW*R0#KrS»N/l60» 1000.0 *32* 15» 1000/ 

7. 

HE**-L)  ( 1 7 ) ( ( 2 ( «J*  I ) » J=1 » 32 ) » 1 = 1 * 128) 

8. 

PI=4.o*ATAM(1.0) 

9. 

P=2*PT 

10. 

LLr'M'U' 

11  . 

UFr 1 . 0/ 1 1704 . 0 

12. 

UO  4 r=l*Ll 

13. 

F ( I ) = ! + ( 1 -LL/2  ) *DF 

14. 

3 

T ( n=r.XP(-(PI*(  1-F(  I ) ) *TW)**2) 

15. 

U=2*Ro*F  ( I ) 

16. 

U=P*(n-INT(U) ) 

17. 

4 

Y(  1 >=t<  I)*S*CLXP(Cr>!PLX(O.0,D)  ) 

18. 

UO  ion  1=1* LL 

19. 

DO  80  J=1*K 

20. 

80 

Z(J* 1 )=Z(J* I)+Y(I) 

21. 

100 

COriTIfiUt 

22. 

WRITE(IB) ( (Z(J*I) »J=1*K  )*I  = 1*  2n  ) 

23. 

STOP 

24. 

END 
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TABLE  3 

THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 

from  copy  furnished  to  DD.Q " 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 

12. 

13. 

14. 

15. 

1 6 • 

17. 

IB. 

19. 

20. 

21. 

22. 

23.  17 

24. 

25. 

26. 

27.  3 

28. 

29.  1 

30. 

31. 

32. 

33. 

34.  2 

35. 

36. 

37. 

38.  5 

39. 

40. 

41.  6 

42.  3' 

43. 

44. 

45. 

46.  4'. 

47.  55 

48.  9* 

49. 

50. 


PROGRAM  Im3 (TAPE 17, T6PE21 i TAPt20. OUTPUT. rAPE6=0UTPUr) 
COMMOM  7 ( K , 2n  HU(lOOO)  tX(lOOO)  ,MC2000) 

OIPCNST^N  F ( 2n  ) » Y ( 2n  ) . SEED ( 3 ) 

DHiEUST  EX  (loo)  *AL  (100) 

COUPLE*  Y , Z , T LMP » TK  * EX 
INJEGEp  S.T* 

Data  tvx  . ro » i*. » s » n/ 

READ ( 1 7 > ( (Z ( J. I ) » J=1 ,32) *1  = 1*128) 

P 1=4.0# A TAN  ( 1 .0 ) 

P=2oPI 

LL=4#K 

Of =1,0/ 1 1704.0 
SEEO(1)=0.0 
SEED ( 2) =1 .0 

CALL  pa*.'3F  ( -20 0 0 * R , SEED ) 

CALL  P»w3F(2000.H|StEO) 

DO  17  t.l.N 
R ( 1 )s»a*1000*R(I>-500 

X(l)a?nnO*R(N*I)-loO0 

AB*H  ( I 1 *»2«-X  ( I ) «#2 
Riiiss^m") 

0 (I ) =X ( T ) /*  (1  ) 

COiMTlN'ir 

CALL  NORMAL (-1000, X,SE0) 

call  noomal(1ooo*x»sed) 
do  3 l-i ,LL 
F(l)ml*fI-LL/2»*0F 
DO  1 1=1.100 

EX(U=rcXP(cUPLX(0.0,-P*(.0l#I-.005)  ) ) 

DO  90  TT=1,LL 
I =LL ♦ l - T I 
DO  30  **=-3Ct32 
DO  2 J=1 . TOO 
AL(J)=n.O 
DO  5 J= 1 * M 

D*F ( 1 ) #(2*R ( J)-KK*0( J) .X (J) ) 

Ms  (D-IMt  f r>>  ) #100*1 
AL  (M)  sAI.  CM)  +1 
Y (KK*31) * C 0.0, 0.0) 

DO  b Js’.lOO 

Y<KK.3M=Y(KK.31)  ♦ AL  ( J)  oEx  ( J) 

CONTlNtir 
DO  5b  KKs)  , K 
TN=(0. 0.0,0) 

DO  40  * 1 a 1 , K 

TK=TK-»7  (Ky  . I ) *Y(KK>M-l  > 

WhlTF  PMTK 
CON  T 1 Ni  ip 
STOP 
END 


r 


— 
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1. 

PKO(>Ka''I  IMlC  (TaPL22,TaP£;23, OUTPUT. TaPi;6  = oUTPUT) 

2. 

CUMMOW  / { K » 2n  ) 10(  1000) * X ( 1 0 0 0 ) ,k(2000) 

3. 

DIMENSION  F<  2n  ) * Y ( ?.n  ) , SLED  ( 3 ) 

4. 

DIMENSION  Ex { 100) »4L (100> 

5. 

COMPLEX  Y.Z.Tf.flP.TK.LX 

6. 

IMTEUEP  S . TW 

7. 

DATA  T*/,Ro»K»$»N/  

8. 

r 

READ ( 2JM  < (Z(J.I) *J=1  *32) *1  = 1*128) 

9. 

PI=4.0#aTAM(1.0) 

10. 

P=2#PI 

11. 

LL=4*< 

12. 

OF=2 • 0 / 1 1704. 0 

13. 

SEED  ( 1 ) =0 . 0 

14. 

SEED ( 2 ) s 1 , 0 

15. 

CALL  pa- 3F (-2000*R*SEED) 

16. 

CALL  RA"3F (2000.R*SEED) 

17. 

DO  17  1 = 1 , N 

18. 

R (I ) =wn*loOO*R (I )-500 

19. 

X(l)=2naO#R(N*l)-l000 

20. 

Atl  = R ( I )•■■>*?♦  x ( I ) «*2 

21. 

R ( I ) = ST'T(Ad) 

22. 

Q(I)=X(1 )/R(I) 

23. 

17 

CUNflNH^ 

24. 

CALL  NORMAL (-100C *X*SED) 

25. 

CALL  NO vMaL ( 1 UOO * a ♦ SED ) 

26  . 

DO  3 I=i .LL 

27. 

3 

F ( 1 ) s[ . ( I -LL/2 ) *DF 

28. 

DO  1 1=1.100 

29. 

1 

EX(I)=c=-XP(cmPLX(0.0.-P<i  (.01*1-. 005))) 

30. 

DO  90  TTsl.LL 

31. 

I =LL+  1 - T I 

32. 

TEMp=(A. O.0.0) 

33. 

DO  2 o kk  = i~*n 

34. 

20 

TEMP=TFMP*Z(^K,I) 

35. 

D=2*K0*r ( I ) 

36. 

D=-P* (d-InT (D) ) 

37. 

TLhP=TF-P*S*CEXP (cMPLX(O.O.d) ) 

38. 

DO  3o  <k-  = -3q,32 

39. 

DO  2 J-i ,100 

40. 

2 

AL(j)=n.0 

41. 

DO  5 J=i«N 

42. 

D=F(I)#(P*R(J)-KK«MJ(J)  ♦X(J)  ) 

43. 

M=(D-INt  <D) > » 1 0 0 ♦ 1 

44. 

5 

AL  (M)  =A(.  (M)  ♦ 1 

45. 

Y(KKOll*(0. 0,0.0) 

46. 

DO  6 J*i *100 

47. 

6 

Y(KK*3l i=y(NK*31)»AL(J)*EX(J) 

48. 

30 

CUNTin,,c- 

49. 

DO  55  kxsi,K 

50. 

TK=(0. 0.0.0) 

51. 

DO  4o  *l*l»i\ 

52. 

40 

TK=TK*7^1#I)*Y(kk*k1-1) 

53. 

TKsT**TFMP 

54. 

55 

WRITE (?^>TK 

55. 

90 

CUNT  Ini ir 

56. 

STOP 

57. 

END  
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TABLE  5 


1. 

PROGRAM  IM61  T AlJE  1 7 » T APE  1 8 . T APE 19  * OUTPUT » TAPE6=0UTPUT ) 

2. 

COV-MOm  2(32*128)  . A(40t>6) 

3. 

Ur'tMslON  Y(  120) *T(128) 

4 . 

DIMENSION  IW(100)  n-.'(lOO) 

5. 

LOGICAL  LE(1()0) 

6. 

equivalence  ( i w ( i ) » w n ) r lf ( i ) ) 

7. 

COMPLEX  Y»Z*TT#TEMP 

8. 

LL=12n 

9. 

UT  J = 1 1 704.0/256.0 

10. 

HE  MM  JM)  ( (Z(J* J ) *Jrl. 32) *1  = 1*128) 

11. 

UO  15  «J-0  * 31 

12. 

DO  5 T=1*12B 

13. 

5 

Y ( 1 ) ~2 ( J+l * 1 ) 

14. 

UO  33  1-1 »LL/2 

15. 

TE  'P-Y  ( 1 ) 

16. 

Y ( 1 ) =Y  t LL/2+ 1 ) 

17. 

33 

Y (l.L/?+l  )=TEMP 

10. 

UO  25  K.K-1 . LL 

4 

14. 

25 

Y ( K K ) rCON JG ( Y ( K K ) ) 

20. 

CALL  Ff-  r P ( Y » LL  * I W * '•/ * LF ) 

» 

2)  . 

UO  37  1=1. LL/2 

22. 

TE’-Pry  ( i ) 

23. 

Y ( D=y(LL/24  I ) 

24. 

37 

Y(LL/?+l)=TEMP 

25. 

UO  10  1=1.128 

26. 

10 

At  128iJ+l)=CAt)5(Y(I)  ) 

27. 

15 

COPT  It.Ut 

20. 

CALL  PLOTS  ( 1 * - 1 * 13HPP0F . MDELl.F.R  * ZQ.  0 * 3HIM2J  . 

24. 

CALL  fKALE  (A  * 4096*  1 * 1 • 0 * AN"IN»  OA  ) 

30. 

UO  45  J=0*31 

31  . 

T(  H:.un*  (LL/2-1  )-UT  1/2 

32. 

UO  18  1=2  * LL 

33. 

IB 

T ( I ) rj ( 1-1 ) +UT1 

34. 

CALL  8KALE  ( T*LL*  1 * 7 • 5 * TMIN »DT2 ) 

35. 

CALL  pLOT(Ttl) *A(120*J+1) *3) 

36. 

UO  20  1=2.123 

37. 

20 

CALL  PL.ort  T(  1 ) . At  1 pfii-J+I  ) .?) 

38. 

CALL  PLOT (0.0* 1.0. -3) 

39. 

4b 

COUTImUE 

40. 

call  plots ( - 1 ) 

41. 

STOP 

• 

42. 

END 

